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Stability Configurations of Electrons on a Sphere 


An unsolved mathematical problem created by J. J. Thomson's efforts to 
Visualize the nucleus of an atom is the determination of stable configurations of 
@ given number (m) of electrons bound to the unit sphere and interacting under 

nutual (newtonian) repulsion. Clearly, one such position would be that of lowest 
potential energy, but uniqueness and presence of symmetries are still unsettled 
uestions. As an analytic procedure a stable configuration can most often be deter- 
nined by a lucky guess supported by the explicitly calculated potential energy, 
thecked to be a semi-definite quadratic form in the infinitesimal displacements 
rom the conjectured configuration. On the other hand from the point of view of 
the computing machine a feasible method is to start with an “arbitrary” con- 
figuration and attain the minimum potential energy by descent. The former 
method (of infinitesimal displacements) was used by L. Féppl [1] to find stable 
configurations for all n < 8 and several larger n, some of which are summarized 
below in terms of “‘rings.”’ 


+1 8=14+3+3+1 
+1 10=1+4+4+4+1 
+1 


12=1+5+5+4+1 
144=1+6+6+1. 


The notation ‘““n = 1 + m + 1” denotes one electron at each pole (0, 0, +1) and 
electrons forming a regular polygon (or ring) on the equator (z = 0). The 
ptation “n = 1+m-+m-+ 1” denotes one electron at each pole and two 
egular polygons of m electrons at equal and opposite latitudes situated so that 

h electron in the upper hemisphere (z > 0) is antipodal to an electron in the 
bwer hemisphere (z < 0). Stable configurations for m = 4, 6, 8, 12, 20 were seen 

iby Féppl to correspond to the five regular solids. The first two omissions in 
Sppl’s work were m = 9 and m = 11. The purpose of this study is to find stable 
tonfigurations for these values of m by descent, using the IBM 701. The principal 

w results (December 1955) are that m = 9 leads to three rings of equilateral 
fiangles with rotational symmetry, while m = 11 is quite irregular, having 

anar but not rotational symmetry, on the basis of the numerical evidence 
ibmitted here. 

We start with an “arbitrary’’ configuration standardized as P;,® = (0, 0, 1), 
P, = (0, .6, —.8), Ps = (.6,0, .8), P?,4 = (sin 2nt/(n — 3), cos 2xt/(n — 3), 0), 
Where 0 < ¢ < n — 4, so as to avoid superimposing any artificial symmetry. 
lore generally let us consider the main loop as consisting of one step in the 

escent, where the general configuration is given by 3m quantities: 


) P; ard (Xi Vis 2), 1 <i _"*% (x? + y? + 37 = 1). 
117 
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The (tangential) force on the particle at P; is (X;, Y;:, Z;) where 





etc., 


eo i — ¥i) — ei — xs)mi + Os — Vidi + (2 — 4))25) 
cieiers = C(xs — xs)? + (vs — ys)? + (@: — 2)? } , 


and the point P; is moved to (x,’, y,’, z:’) where 
(3) xi! = x; + AX;, etc. 


Here h is a positive constant to be determined automatically as described later on. 
Then, replacing P; by Pi(h) = (x;(h), y:(h), 2:(2)), where 


(4) xi(h) = x;'/(xi/2 + y/? + 2), etc., 
we close the loop. The exit occurs when 
(S) € = max {|x; — xi(A)|, |¥s — vi(A)|, |2s — 2i(h)|} <e. 


The choice of “‘step-size’’ h has to be made internally. We consider only those 

h of the type 4 = hop* where hp is a positive parameter (initially .1 and variable 

with each circuit of the main loop), p(>1) is a fixed ratio taken as 2!, and is 

an integer to be chosen. We call V(= > |P:P;|-") the potential of the configura- 
i>j 


tion P; and V(h) that of the configuration P;(h). There are two possibilities. The 
first possibility is that V(’o) < V. In that case we choose hk = hy (k = 0) urtless 
V(hop) < V(ho), in which case we choose & to be the minimum positive integer 
for which V(hop') < V(hop*), ¢ = 1, 2, ---, &, while V(Aop*t") > V(hop*). (We 
must exclude (by programming) the possibility that # be © or that the 3n 
corrected coordinates x;,’, y;’, 2; be parallel to the force vectors X;, Y:, Z;. This 
might happen accidentally when radius vectors parallel to the force vectors repre- 
sent a configuration of lower potential energy than the configuration P;.) Then 
h = hop* and h replaces fo. The second possibility is that V(’o) > V. In this 
case, letting & be the minimum positive integer for which V(hop-*) < V, then 
h = hop with h now replacing ho. (Here again 4 might become zero through 
round-off loss in the computation of V(h), and this possibility must be excluded 
by programming.) 

In practice, h = hp in at least 3 of the cases; and when / charged, it changed 
only by the factor p*!. Generally .2 <h < .4. 

The program was coded in Speed Code III and read in from IBM instruction 
cards and floating decimal cards for n, h, p, «, etc. The use of a binary instruction 
deck cut the read-in time of 150 cards per minute to less than two minutes. The 
machine generated the P; internally and, in cases of reruns, accepted the last 
P; and h and modified entrance instruction from correction cards. Each major 
(descent) loop took about .3n? seconds (slightly longer if k was changed by the 
loop). At the end of each loop the machine loaded the pairs h, V(h) in an output 
block of the memory, dumping the pairs 25 at a time to provide some means of 
observation. The machine also loaded the 3n + 3 values ¢, V, h, P; on tape, 
although there was no occasion to dump the tape later on. As a further monitoring 
device the sensing switch P was made available to print these 3n + 3 values at 
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the end of each loop when desired. At the exit, the residual (h, V(#)) pairs, 3n + 3 
aforementioned values, and the m* mutual distances |P;P;| were printed out; 
and the Q switch was used to bring the computation to the exit at the end of the 
current loop in case of shortage of time. Print-out time was used sparingly since 
each printing operation takes 5/3 + 21/5 seconds for / lines of 5 words per line. 

The descent described becomes badly oscillatory; and although we need not 
regard this fact as wholly unfavorable, we should see the descent in more conven- 
tional terms. (The oscillations simulate vibrations in that normal modes can be 
identified. This will be the subject of a more complete study later on.) We 
imagine the particle positions parametrized in terms of 2m independent variables 
qi such as two of the three displacement coordinates of each point from equi- 
librium. We assume small displacements. Then some real, symmetric, and (pre- 
sumably) positive, semi-definite matrix ||a;;\| exists such that 2V = ¥ aiqaq;. 
The matrix has the eigenvalues 4;(>0) and normal coordinates Q; for which 
2V = > A.Q2. Thus equation (3) becomes approximately 


(3a) Aqi = —hdV/dqi = —h>d aisgi, 
or 
(3b) AQ; = —hQi, 


suggesting exponential decay rather than ‘“dynamic’’ oscillation. The oscillatory 
behavior here is (presumably) caused by the Jarge range of non-vanishing A;, or 
by the fact that no one hk can make ail the ‘decay ratios” | (Q; + AQ)/Q:| 
=|1 — hAd;| appreciably less than unity (for A; > 0). 

Certainly, for purposes of establishing a (conjectured) local minimum for V, 
it would be reasonable to check the A;. We are concerned more with the pre- 
liminary process of obtaining the initial conjecture in the face of very slow con- 
vergence. The procedure will be to cut the number of degrees of freedom, and 
(effectively) the range of A; through conjectured symmetries. For instance when 
n = 9 (m = 11) for the initial configuration P; described earlier V = 27.07665314 
(= 43.03629491) and after 32(51) circuits of the main loop V had descended to 
the value 25.76006543 (40.59907213), stable to four (two) decimal places, while 
the oscillations in particle position, ~, had the order of magnitude of 10-* (10-*) 
with no perceptible improvement. Yet a drawing on graph paper of the “‘final,” 
ie., the thirty-second (fifty-first) approximation for m = 9 (m = 11) displayed 
enough symmetry to enable us to cut down the number of degrees of freedom 
from 18 to 1 (22 to 5). 

In particular, when m = 9 it becomes graphically clear that the arrangement 
is9 = 3 + 3 + 3, ie., an equilateral triangle on the equator with two symmetric 
equilateral triangles, one in each hemisphere, giving the nine points (1, 0, 0), 
(—43, +3#/2, 0), ((1 — a*)#/2, + (3[1 — a*))#/2, +a), (—[1 — a*}', 0, +a), where 
a is the remaining degree of freedom. 

Likewise, when m = 11 it becomes graphically clear that the arrangement is 
five points on the equator and three in each of the hemispheres, symmetric with 
respect to the equator (z = 0) and meridian (y = 0), yielding the configuration 
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(—1, 0,0), (—a1, + (1—a;’)!, 0), (a2, + (1 —az*)!, 0), (—as, au, + (1—a3*—a,’))), 
((1—as?), 0, +a), with five degrees of freedom a;. 

Thus we observe the final configuration for m = 9 (m = 11) and regenerate tak- 
ing initially a = .699929 (a; = .505694, a2 = .507286, a; = .089001, a, = 492006, 
as = .571672). These are crude guesses based on averages of diagonal lengths. 
In fact V now increases to 25.76028061 (40.64744085) which brings us back to 
the twenty-third (eleventh) step of the descent process! We are further ahead, 
however, by virtue of the new self-preserving symmetry. In fact in five (thirty) 
additional circuits of the main loop we find V has come to the even lower value 
25.75998651 (40.59645048) probably correct to eight decimal places at which 
point the particle positions change by — < 10~-* = e. The final values of the coordi- 
nates of the configuration are given by a = .703648 (a1 = .515358, az = .552626, 
as = .168322, ag = .488923, as = .591930), probably correct to five decimal 
places, in a total of about 15(60) minutes of computing time. 

The author is indebted to G. Pélya and J. L. Ullman for illuminating discus- 
sions and to Don Hart and George Ryckman of the General Motors Technical 
Center at Warren, Michigan for kindly making the IBM 701 available. 
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An Iterative Method for Taylor Expansion of 
Rational Functions, and Applications 


1. The iteration and some applications. A simple iterative procedure can be 
set up for determining the coefficients of the Taylor expansion of a rational func- 
tion at any point in the complex plane other than one of the singularities of the 
function. This procedure would have many applications, particularly in the 
evaluation of inverse Laplace transforms. We shall discuss the general method 
first, and then some applications, followed by a discussion of truncation and 
round-off errors. 

Suppose f(s) is a rational function, that is, 


- Ps) 
(1) H6) = 


where p(s) and q(s) are polynomials. We shall assume that we have p(s) and 
q(s) expressed in powers of s, rather than in their factored form. 
Now consider any point ‘‘a’”’ in the complex plane, with the restriction that 


q(a) # 0. Then by Taylor’s theorem, within the circle of convergence around “a”, 
we have 


rr rong a Rp 
(2) f(s) q(s) 0 + ax(s — a) + a2(s — a)? + ---. 
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Multiplying (2) by q(s), 
(3) p(s) = g(s)[ao + ai(s — a) + a2(s — a)? +--+]. 
Setting s = a, we see that 


(a) ay = 2E), 


Having thus determined ao, we define the polynomial 


(5) 1(s) = p(s) — ang(s). 


Clearly ¢:(a) = 0, so “a” is a root of ¢:(s). Thus we can define another poly- 
nomial, 


_ ofS) 
(6) p(s) = 
Then from (3) and (5), we have 
(7) pi(s) = q(s)[ar + a2(s — a) + as(s — a)? + ---]. 


Comparing equation (7) with (3), we see that we have replaced p(s) by :(s), 
and a; by aj41 (j = 0,1, 2, 3, ---). Continuing as before we get 








_ ?1(@) 

(8) a) q(a) ’ 
(9) $2(s) = pi(s) — aig(s), 
(10) pr(s) = 2) , 

s-—a 

_ p2(a) 

(11) r q(a) , 
and so on 


More precisely, defining po(s) = p(s), we get the recursive relationships 





p;(a) 
1 a; = ’ 
- ‘= "q@) 
(13) $i41(S) = p;(s) — asg(s), 
and 
(14) . Sinale) = 2, 


s a 


where j = 0, 1, 2, 3, etc. 
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This method is particularly suitable for digital computation for several 
reasons. One is that since the process is recursive, it is just as easy to program 
for a large number of the a; as for only a few. Also let us examine the equations 
(12), (13) and (14) in more detail. We start with the coefficients of the poly- 
nomials o(s) (i.e., p(s)), and g(s). We can obtain fo(a) by performing the syn- 
thetic division of po(s) by (s — a) and taking the remainder. We obtain g(a) by 
the same method, and then perform the straight division to get ao. The coeffi- 
cients of ¢1(s) are easily calculated, since they are obtained by subtracting ay 
times the coefficients of g(s) from the corresponding coefficients of po(s). We now 
again use our synthetic division process, dividing ¢:(s) by (s — a), taking the 
quotient for :(s) and inserting an automatic check to make sure that the re- 
mainder is zero (to within round-off errors, etc.). Then we go back to equation 
(12), and compute a: = p:(a)/gq(a). We evaluate :(a) by the synthetic division 
method again and of course g(a) is now already calculated from the first round, 
and does not have to be calculated again. Thus the calculation continues, and a 
counter can be used to control the number of iterations, depending on the number 
of coefficients a; that is desired. 

Thus we have reduced the problem primarily to synthetic division. Of course 
in general we will have to use sub-routines for complex arithmetic, although if 
the original polynomials have real coefficients and we are only interested in 
expanding about points on the real axis, this will not be necessary. 

We now consider some applications of the above method. First suppose we 
wish to find the inverse Laplace transform of a rational function 


(15) P(s) = wat 





where N(s) and D(s) are polynomials, the degree of D(s) being higher than that 
of N(s). Suppose also that we know the roots of D(s), so that 


D(s) = K(s — s1)™(s — S2)"* +++ (S — Sm)™. 


Then we can expand F(s) in partial fractions: 








Au Ais Ain, An 
(16) F(s) = peas + rage ae boa 
+ _ Ar fee _ Armin _ 
(s — Se)? (s — Sm)™™" 
or in contracted form 
m nm Ay 
(17) F(s) = XL 


imi jor (S — Si)? 
Once we have F(s) expanded in this fashion we can easily find the inverse trans- 


1 tle 

j — A. » 5 3 ++. . 
[im @-in” ) 
(0! = 1). The problem is to find the coefficients A;;. Clearly if we can find the 
coefficients corresponding to 51, i.e., A1;, we can find the rest in the same manner. 
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Thus let us break up the expression of equation (17) into the contribution from 
the root s1, and the remainder. We get 


(18) P(s) = & + Ax erie Ain, 


—s; (s—s;)? (s — s1)™ 





+ h(s), 


and h(s) is analytic at the point s;. Returning to equation (15) now, we write 


N(s) 


PO) = Ga a)Die)’ 





where the polynomial 


— 
(20) D,(s) it (s <a $,)™ 
Clearly D,(s1) # 0. Then multiplying (18) by (s — s1)™, and using (19), we get 


N(s) 
D,(s) 





(21) = Ain, + Axgn,—y(S — $1) + Arqn,—m(S — 51)? + --- 


_ Ai(s - $;)""—* an Aun(s - $3)" + (s - $;)™h(s). 


Thus Ain,, A1gn,—»), *** Aiz, Ans, are the first ; Taylor coefficients for the expan- 
N(s) 
D 1(s) 
described. We would wish to have D;(s) in powers of s rather than in the factored 
form, so we would start with D(s) expressed in powers of s. Then we perform the 
synthetic division by (s — si) m; times to get the coefficients of D,(s). The 
remainder should be zero each time, and as before this fact can be used for 
error checking. 

In most cases, F(s) will have only simple poles; i.e., we would have m; = m2 --- 
= mm = 1. In that case the method described by Titus [1] is essentially the same 
as this one, perhaps even a little simpler. However, Titus’ method gets consider- 
ably more cumbersome when treating multiple poles, and since the method we 
have described here works equally well for multiple poles, it would seem prefer- 
able in this respect, provided the computing machine used has sufficient internal 
memory for the incorporation of the necessary subroutines. Furthermore the 
routine for the Taylor expansion of rational functions could be used for other 
applications such as will be described below, as well as for partial fraction ex- 
pansion. 

In the next example we again consider 





sion of around the point s:, and we can find them by the method we have 


(15) P(s) = ae 


where the degree of N(s) is lower than that of D(s). However, suppose we wish 
to have the inverse Laplace transform, F(#), as a power series in t. To do this we 
expand F(s) around the point at infinity, using a method due to Heaviside [2]. 
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First we substitute s = 1/z, and then multiply numerator and denominator by 
2", where N is the degree of D(s), obtaining 


1) _ «) . 2@ 
(22) p(+) == 0 


where p(z) and g(z) are polynomials. Also (due to the fact that N(s) is of lower 
degree than D(s)) 


(23) f() =9, 

so that we can factor out z and write 

(24) f(2) = 2g(2). 
Then if 

(25) e() = E ap, 

(26) F(t) = =e (t > 0). 


Thus we can use the method we have described to obtain the coefficients a; in 
(25), and substitution in (26) will give the series for F(t). The rate of convergence 
of this series will be discussed later in some detail under the subject of trunca- 
tion errors. 

The formulae for a; in this case can be simplified somewhat. Suppose the 
function F(s) is in the form 


aosN—! + aysN—-2 4+ .-- + ay_os + aw_1 4 
SN + bisY—1 + bos¥—? + --- + by_is + by 


(It can always be expressed in the above form, with the coefficient of the highest 
degree term, s¥, equal to unity.) Then from (22) and (24) we obtain 


(27) F(s) = 





(28) g@) - 

where 

(29) Po(2) = ao + aye + «++ + ayy", 
and 

(30) Q(z) = 1+ bis + bee? + --- + bys. 


It is not difficult to show from (12), (13), and (14) that the coefficients a; in 
(25) are given by 


(31) aj = 0,3; 
(32) Op, i41 = Oe415 — Odinr, OLR N—2, 


(33) Gn-1,j41 = — ajby, 
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where 
(34) Ano = Gm, OS mS N—1. 


One advantage of the above method is that it is not necessary to know the 
roots of the denominator of F(s). 

The third application we shall mention is one from probability theory. Suppose 
we have a random variable X, which can take on the values 0, 1, 2, 3, etc. (and 
no other values) and that it takes on these values with probability Po, P:, P2, P:, 
etc. Thus P; (j = 0, 1, 2, ---) gives the “probability distribution” of the random 
variable X. 

In dealing with such distributions, especially when treating more than one 
of them at a time, it often is convenient to treat not the probability distribution 
itself, but its “generating function,’ P(s). The generating function is precisely 
that function which, when expanded in its Taylor series around s = 0, has Taylor 
coefficients Po, P:, P2, etc. Furthermore, these generating functions are often, 
though not always, rational functions. In such cases we can see that the method 
we have described would be very useful in obtaining a probability distribution 
from its generating function. 

A good discussion of generating functions is given in [3], chapter 11. 

There are also many other possible applications for the procedure that has 
been presented here. For example it might simply be used to calculate one of the 
derivatives of a rational function, since the m-th derivative at a point “a” is 
given by 


(35) EE I. ik 





where a, is as before coefficient of (s — a)" in the Taylor expansior about “‘a’’. 

2. Truncation and round-off errors. We come next to the important question 
of errors in the above method. We shall consider first truncation errors, due to 
terminating the series (2) after a finite number of terms. To do this we assume 
that each term in the series is evaluated exactly, without round-off error, leaving 
the discussion of round-off errors till later. 

First we note an important theorem. Let p be the radius of convergence of 
the series (2) about “‘a.’’ Let 0 < ro < p, 0 < 6 < 1, and let r; = Oro. Let M be 
an upper bound of | f(s)| on the circle, |s—a| = ro. Then if the remainder 


(36) Ra(s) = E als — a, 
Me 
(37) |Ra(s)| < er 


whenever |s — a| < 1. For a derivation of this condition see [4], chapter VI. 
This gives uniform convergence on the region, |s — a| < 11, and we can 
make | R,(s)| less than any given positive number y, over the region, by choosing 
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n such that 





(38) n> 


where ‘“‘In’’ denotes the natural logarithm. 

In the first example we have discussed, pertaining to partial fraction expan- 
sion, the question of truncation errors does not arise, since in this case we are 
interested only in a finite number of the terms. In the second example however, 
the power series expansion for F(#), the truncation errors are important. Here we 
cannot use (38) directly because we are interested in the truncation errors in 
equation (26), rather than (25). In order to study this case we first define 





(39) h(z) = g(z) — g(0) 
z 
Then 
a P,(2) 
(40) h(2) “a Q(z) ’ 


where Q(z) is given by equation (30) and 


(41) Pi(2) = Co + crt + +++ + Cy_s2", 
where 

Co = G1 — dob 

C1 =de- Aobe 2 
(42) 


Cn-2 = Gn-1 — Gobn-1 
Cv-1 = — doby. 


(The coefficients c, above are the same as the quantities a;,1 in equations (32) 
and (33).) The MacLaurin series for 4(z) is given by 


(43) h(z) = X By" 
i-0 
and from (25) and (39) we see that 
(44) B; = Aj+1 Gj 9 0, i. 2, pe. -). 


Next we let 7: be a positive number, as yet unspecified except that it is less 
than the modulus of each of the zeros of Q(z). Then (43) is uniformly convergent 
for |z| <r. 

Now suppose we wish to take enough terms of (26) to be accurate within +e 
on the interval 0 < ¢ < t%. We can find mp such that when n > mo 
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for |z| < 71, where 


(46) Ra(s) = 5 Be. 
j=n 

Then, if 

(7 Fa) = 5, 


we have, by the inversion theorem, using (22), (24), (25), and (44), (¢ > 0), 


sit 


ws) 1FO-FOl= Z| Poas-E af Sas 
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nv~ ern, (7) (+4) Ji 
on f eR, P a Te idy 
. 7 


~ a - = (445) 
T1 7, 


1 2 io" 2 dy 
(49) |F@® — F,()| +s * :. i-* 


ds 











wae s? 
"1 


Thus if 0 < ¢ < ft, and ” > mo, we obtain 





The above approach is equivalent to that used in [2], chapter XIII to 
prove the validity of equation (26), in the more general case where F(s) is 
analytic at infinity, but is not necessarily a rational function. 

We now wish to make a specific choice of r; and an appropriate mp». In order 
to use the theorem stated above for truncation errors, we first pick ro and then 
make 7; < ro (both being less then the radius of convergence of (43)). First let B 








be the maximum of the quantities |b:|, |b2|, ---|bv|, obtained from equation 
(30). Then we will choose 
1 
(50) ro-TaB 
For |z| < ro, 
Seer aa 
= Ss i+B\i+B 1+B 
N 
+ trdos) 
_— B ee 1 
+P, _/( ..) (+ By" 
1+B 


so that our choice of ro is a satisfactory one. 
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Next we must find M, an upper bound, of |h(z)| on the circle |z| = rp. 
From (40), (41), and (42), we see that each of the coefficients in the numerator, 
P,(2), is less than A(i + B), where A is the maximum of the quantities |a;|, 








1 
|a2], --+|@w-1|, so that for |z| = ro = vet 

2) |P <aats)|1 1 we (5 
(52) [Px@)| $40 +B) 1+ 5+ (p55) +: +l | 
mrs J 
= A(1 + B) 

pale: 

1+B 


Then using (51), we have for |z| = ro 
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1 EE 
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(53) |h(s)| < AQ + B)# 


so we choose the right hand side of (53) for M. Letting 


6 
(54) nen o <6 < i, 
(45) becomes 
(1 B _(1+Byto 
(55) ae <2. o, 


6 


Thus, from (38), (55) will be satisfied for |z| < 71 if m > mo, where 
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and A and B are the maximum absolute values of the quantities a; and b; respec- 
tively in equation (27) (WN being the degree of the denominator of F(s)). 
Thus if m > mo, defined by (56), we will have 


|F() — Fa)|<« O <t# <b). 
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6 can be any number between 0 and 1, but we might choose it so as to make 
the term in ¢ a minimum. If we do this we obtain @ = ¢~', so that (56) becomes 


4 + BUG +BY ~ 13). 
2B(e — 1)e 





(57) No = e(i + B)to + in| 
Finally we consider the third example of the probability generating function 


(58) P(s) = ¥ Ps? 
7=0 


Since we must have >> P; = 1, the radius of convergence, p, is either greater 
j=0 


J 
than or equal to unity. The probability that X > n is 
(59) Ra(1) = L Pj. 
j=n 


One might often be interested in choosing m large enough so that this quantity 
is less than some preassigned small number. If p > 1 this can be done by the 
general theory we have described. If p = 1 however we see that this theory does 
not apply and the problem becomes somewhat more difficult. a 

We now come to the general problem of numerical round-off errors. If one is 
summing a certain number of terms of a series and one can estimate the accuracy 
of each term, it is not difficult. to study the error in the sum. The more serious 
problem, in our case, appears to be the possibility of cumulative errors in the 
calculation of the a; themselves. 

Let us assume that we are using a floating decimal computing system which 
retains four significant figures in the mantissa. Consider the calculation 


(60) x = [4.790 — 4.751] + 1.783 


where we perform the subtraction first, and then the division. The computer 
will give 


(61) x = 2.187 X 107. 
Next consider the equivalent calculation 
(62) x’ = (2.507 XK 4.790 — 2.507 X 4.751) + (2.507 X 1.783) 


where we evaluate each term in the numerator, subtract, evaluate the denomi- 
nator, and finally divide. The computer will give 


(63) x! = 2.237 X 107. 


In this case we could see beforehand that the answer would only be valid to two 
significant figures, and we note that x and x’ agree to two figures but no further. 

This illustrates the point that round-off errors can usually be detected by a 
duplicate calculation by a modified routine, although it is difficult to obtain hard 
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and fast rules for doing so. In the above example, for instance, if we substitute 
some other number, k say, for 2.507, the desired effect does not occur for all 
values of k. An integral power of ten will not work (since we are assuming round- 
ing-off is done in the decimal system). Also one finds experimentally in this case 
that k = 1.001 and & = 2 are unsatisfactory choices. 

In our present study we might use a transformation of the type 


(64) s = ks’ 
and re-evaluate the a; in equation (2) by expanding the function 
(65) o(s’) = f(ks’) 


around the point 


ia? 
(66) a a 


obtaining quantities a,’ as coefficients of (s’ — a’), from which, theoretically 


at a 
(67) a =F; 


Then by comparing the original quantities a;, with those obtained from (67), as 
j increases, one could estimate the cumulative building up of round-off errors. 
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Numerical Integration over Simplexes and Cones 


1. Introduction. In this paper we develop numerical integration formulas for 
simplexes and cones in m-space for m > 2. While several papers have been written 
on numerical integration in higher spaces, most of these have dealt with hyper- 
rectangular regions. For certain exceptions see [3]. Hammer and Wymore [1] 
have given a first general type theory designed through systematic use of cartesian 
product regions and affine transformations to extend the possible usefulness of 
formulas for each region. 
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Particular formulas were developed by Hammer and Wymore for certain 
symmetrical type regions including spheres and hyperspheres, cubes and hyper- 
cubes. The methods they have used make it possible to obtain numerical integra- 
tion formulas for a much larger class of regions than heretofore. However, it is 
not practical to obtain specific integration formulas for all regions of interest. 
Hence, in this paper we develop methods for obtaining numerical integration 
formulas over simplexes which may be used, in principle, to approximate other 
regions. Since one of the methods is applicable to cones in general if a formula is 
given for the base, we include that in the development. 

While we have calculated certain formulas in specific cases for n = 2, 3, and 4, 
we use these as illustrations and hope to build a more complete table later. 
Inductively we give integration formulas holding exactly for the k-th degree 
polynomial in m variables over the n-simplex. Another method which may be 
more valuable in some applications requires affine symmetry of the evaluation 
points. Here the general theory is missing, but a general type of method is pro- 
posed for which illustrations are provided. 

2. Approximate integration formulas for cones. In the following a theorem of 
Hammer and Wymore will be used which we state as follows: 

THEOREM 1. If 


Ease) — f sav = BO, 
then 
E Wagln) — f eindav = WEL, 


where T is an affine transformation, » = AE + no, of E, onto itself; g(n) = f(£); 
W is the absolute value of the determinant of A; R is an n-dimensional region included 
in the domain of f and &1, ---, &, «++, are points in the domain of f. 

This theorem permits us to consider convenient, specific regions R to develop 
formulas for the class of all their affine transforms. 

We want integration formulas of the form 


(1) f fav = Daft), 


where the numbers a; are constants; £1, ---, & are points in the domain of f; R is 
a bounded closure of an open set in Z,. From this standpoint the most interesting 
simple regions are the simplexes, which are special cones, since every polyhedron 
is composed of a finite number of simplexes. 

Let an n-dimensional region R be embedded in the hyperplane x = 1 in E,,; 
where we represent the points in Z,,; by (&, x), where é is a point in E,. Then the 
set of all points xR, where 0 < x < 1, is a cone C with base R and vertex at the 
origin in Eq41. 

Let f(é, x) be a function defined over C and suppose that a suitable numerical 
integration formula is given over the base R of C. If, for example, 


Q) [se nav. =D asee 0, 
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then 
@ fi seqnav= fae f senav.= f x E asflats 2d, 


since the Jacobian of the affine transformation from R to xR is x—*. Define a 
function 


g(x) = zu ajf (xé;, x) 


f fdv = i ; xg (x) dx. 


Now we ask for numerical integration formulas of the form 


and then we have 


(4) f * eng (x)dx = X beg(x:). 
0 i 


Since such formulas may certainly be found we then have 


(5) Jf sa0 = 2X bases 20, 


which is of the form required, that is, a weighted sum of integrand values. 
THEOREM 2. If a formula (2) holds precisely for polynomials f of n variables of 
at most degree m over a region R and if a formula (4) holds precisely for polynomials 
g of at most degree m in x, then (5) holds over C for all polynomial functions f contain- 
ing terms of at most degree m in its n + 1 variables. n 
We will forego a formal proof of Theorem 2 since it follows from the affine 
invariance of the class of polynomials of at most a certain degree, so that 


x” asf (ety x) = f IG 2)dV 


In order to develop specific formulas one may start with any of the numerous 
formulas for a line segment integration and proceed to obtain a formula for a 
triangle; using this obtain a formula for a tetrahedron and so on. This is the type 
we consider now. In higher spaces, to keep the number of points smaller, the 
integration formulas we choose for /o! x"g(x)dx are based on orthogonal poly- 
nomials with weight function x" so that with m values of x; the formula holds 
precisely for a polynomial of degree at most 2m — 1 by formula (3). For n = 0 
this gives the Gauss integration formula on the line with points of evaluation the 
roots of the m-th Legendre polynomial. Then if the formula (2) for the base is 
exact for polynomials of at most degree 2m — 1 in the first variables, the final 
formula holds for polynomials of degree 2m — 1 at most. Thus by this means 
we obtain with m”* points a formula valid for all (2m — 1)-degree polynomials over 
a simplex in m-space. While we have discussed a range from 0 to 1, the affine 
invariance permits us to give the formulas as valid for all simplexes, the necessary 
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adjustment of coefficients and points being given by the transformation as in 
Theorem 1. 

The orthogonal polynomials, g,,(x), over the interval (0, 1) with weight func- 
tions x" are given by Christoffel’s formula in terms of the Legendre polynomials, 
P,,(x), on the interval (0, 1): 


P(x) Pings (x) + Pan (x) 

P,,(0) Pn+i(0) * Pin(0) 
(6) x"m(x) = | P’n (0) P'm+i(0) “m+n (0) 

PE-” 0) Dacre 0) - PSz? (0) 


In this formula and in the tables the value of m is one less than the dimension 
of the space—for the line, » = 0; for the plane, m = 1; and so on. The resulting 


nNo— 


w — wm 4 w 
wn whe -— ar Whe whe one 


Ch eel 


x; 
0.66666 66666 66666 667 


0.35505 10257 21682 190 
0.84494 89742 78317 810 


0.21234 05382 39152 944 
0.59053 31355 59265 287 
0.91141 20404 87296 044 


0.13975 98643 43780 552 
0.41640 95676 31083 175 
0.72315 69863 61876 278 
0.94289 58038 85482 299 


0.09853 50857 98826 4273 
0.30453 57266 46363 885 
0.56202 51897 52613 862 
0.80198 65821 26391 897 
0.96019 01429 48531 218 


xj 
0.75000 00000 00000 000 


0.45584 81559 88774 713 
0.87748 51773 44558 620 


0.29499 77901 11501 618 
0.65299 62339 61648 121 
0.92700 59759 26850 269 


0.20414 85821 03227 136 
0.48295 27048 95632 480 
0.76139 92624 48137 593 
0.95149 94505 53002 709 


TABLE I. n = 1 
bj 
0.50000 00000 00000 000 


0.18195 86182 56022 831 
0.31804 13817 43977 169 


0.06982 69799 01454 1224 
0.22924 11063 59586 248 
0.20093 19137 38959 648 


0.03118 09709 50008 0822 
0.12984 75476 08232 439 
0.20346 45680 10271 322 
0.13550 69134 31488 149 


0.01574 79145 21692 2766 


Gm (x) 
V4(2 — 3x) 


V¥6(3 — 12x + 10x*) 


V8(4 — 30x + 60x* — 35x%) 


¥10(5 — 60x + 210x* — 280x* 


+ 126x*) 


0.07390 88700 72616 6584 ~12(6 — 105x + 560x* — 1260x* 


0.14638 69870 84669 768 
0.16717 46380 94369 395 
0.09678 15902 26651 7818 


TABLE II. 2 = 2 
bj 
0.33333 33333 33333 333 


0.10078 58820 79825 431 
0.23254 74512 53507 905 


0.02995 07030 08580 6981 
0.14624 62692 59866 020 
0.15713 63610 64886 615 


0.01035 22407 49918 0652 
0.06863 38871 72923 0663 
0.14345 87897 99214 191 
0.11088 84156 11277 894 


+ 1260x* — 462x5) 


Gm (x) 
V5(4x—3) 


V7 (15x? —20x+-6) 


V9 (56x? — 105x*+-60x — 10) 


V11 (210x*—504x3 


+-420x*— 140x-+15) 





0.14894 57870 52983 580 0.00411 38252 03099 00782 
0.36566 65273 69113 217 0.03205 56007 22961 9169 
0.61011 36129 34480 701 0.08920 01612 21590 0168 
0.82651 96792 28304 566 0.12619 89618 99911 440 
0.96542 10600 81784 870 0.08176 47842 85770 9715 


V13(792x' —2310x4 
+2520x* — 1260x*-+-280x —21) 


wn 
ar whe 
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TABLE III. 2 = 3 


j Xj b; 9m (x) 
1 0.80000 00000 00000 000 0.25000 00000 00000 000 V6(4—5x) 

2 1 0.52985 79358 94884910 0.06690 52498 06888 7467 V8(10 —30x+21x*) 
2 0.89871 34926 76543 662 0.18309 47501 93111 252 
1 0.36326 46302 16511947 0.01647 90592 82671 7230 

3 2 0.69881 12691 63613 535 0.10459 98975 56806 681 V10(20 —105x+168x*— 84x?) 
3 0.93792 41006 19874 523 0.12892 10431 60521 608 
1 0.26147 77888 30889 686 0.00465 83670 60069 48897 

4 2 0.53584 64460 88250 229 0.04254 17241 42766 6674 V12 (35 —280x-+-756x2—840x* 
3 0.79028 32299 69286 800 0.10900 43689 38641 000 +330x*) 
4 0.95784 70805 66118 662 0.09379 55398 58522 9295 


TABLE IV. Integration Formula over Triangle Exact for 7th Degree Polynomial 


Ve a 


0.86113 63115 94052 580 0.34785 48451 37453 860 
0.33998 10435 84856 264 0.65214 51548 62546 143 
— 0.33998 10435 84856 264 0.65214 51548 62546 143 
— 0.86113 63115 94052 580 0.34785 48451 37453 860 


PWN &w 


Points of evaluation (x;, x;yx) 


j Xj XiY1 = —XjVa XjV2 = —XjY3 
1 0.13975 98643 43780 552 0.12035 22940 89888 328 0.04751 57045 30876 4547 
2 0.41640 95676 31083 175 0.35858 53991 82305 152 0.14157 13593 61934 441 
3 0.72315 69863 61876 278 0.62273 67399 39136 723 0.24585 96668 98990 366 
4 0.94289 58038 85482 299 0.81196 18147 75453 378 0.32056 66993 96768 242 
Wy = ba, = weights for points (x;, x;yx) 

j k=1,4 k.= 2,3 

1 0.01084 64518 21050 5090 0.02033 45191 28957 5733 

2 0.04516 80985 64739 8624 0.08467 94490 43492 5770 

3 0.07077 61357 96171 8794 0.13268 84322 14099 443 

4 0.04713 67363 86764 6765 0.08837 01770 44723 4729 


Qm(x) are orthogonal, but not necessarily normal. The roots are the values of x; 
required. Normalization of the gm(x) gives the weights 5; by: 


(7) bjt = x Los (xs) P 


where x; is a zero of gm(x). 

The orthogonal polynomials, g(x), over the interval (0, 1) with weight func- 
tion x" are the Jacobi polynomials under the linear transformation x’ = $(1 + x). 
The standard definition of the Jacobi polynomials gives the weight function as 
(1 — x)@(1 + x)* and the interval of orthogonality as (—1,1). In this case, 
a = 0,8 = mand the linear transformation given above takes (1 + x)" into (2x’)* 
and (—1,1) onto (0,1). From this fact, explicit representations of g,,(x) are 
available. 

We consider an example of how to combine a formula for » = 0 (the Gauss 
formula on the line) and one for m = 1, to obtain a formula for the triangle. In 
the following discussion we consider a plane triangle with vertices (0,0), (1, 1) 
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and (1, —1). Use of the formulas may be made for arbitrary triangles by applying 
theorem 1. If we take, in each of the above cases, m = 4, a formula is obtained 
which is exact for the general polynomial of degree 2m — 1 = 7. For m = 1, x; 
shall denote the roots of g.(x), and b; the corresponding weight. The roots of the 
Legendre polynomial, P,, shall be denoted by »,, and the corresponding weights 
by a,. Table IV gives the values of x;, b;, and yz, a,; the points, (x;, x7), at which 
the integrand is to be evaluated, and the weights, wy, at these points. 

The first three tables give: 1) the polynomials g,,(x) obtained from equation 
(6) for m = 1, 2, 3 respectively, and values of m indicated in each table, 2) the 
roots x; of g(x), and 3) the values of b; obtained from equation (7). The calcula- 
tions were made at the University of Wisconsin Numerical Analysis Laboratory 
using a CPC with Eugene Albright’s 18-digit floating decimal board. The approxi- 
mate error in the x,’s is no more than 1 in the seventeenth significant figure; the 
approximate error in the 0;’s is no more than 1 in the sixteenth significant figure. 
The sum of the 3,’s is 1/(m + 1) where 1 is defined as in equation (6). 

3. Symmetrical formulas. While the foregoing development makes it possible, 
in principle, to obtain numerical integration formulas for the m-simplex to hold 
exactly for polynomials of at most degree k, we do not mean to suggest that such 
formulas are the most desirable. One feature of these formulas is that they are 
unsymmetrical—i.e., in a given simplex the particular points of evaluation are 
not an invariant set under affine transformations taking the simplex onto itself. 

In this section we give the preliminary results of our investigations resulting 
from a requirement of affine symmetry. We give these results since we believe that 
formulas for triangle and tetrahedron will be most useful and we have certain 
specific formulas for these regions. 

The requirement of affine symmetry we make is simply this: If an integration 
formula involves calculation of the integrand at a certain point P to be multiplied 
by a weight, w, then all images of P under all affine transformations of the region 
onto itself will appear in the formula with the same weight, w. While such a re- 
quirement would appear to increase the number of points, in dealing with poly- 
nomials this is not the case; actually we have obtained fewer points than with 
the other method. 

Let us represent points in the space as vectors and write the vertices of a 


3 
triangle as Vi, V2, V3, and the centroid as C = 4 >> V;. The first affine invariant 
1 


formula for the triangle is to use the centroid as the sole evaluation point with 
weight equal to the area. This method works for all bounded regions in all finite 
dimensional spaces to give a formula for the general linear function over the region. 
For simplexes or for hypersquares, this is likely to be useful. 

The quadratic function over the triangle we have shown to be integrated 
exactly by evaluations at rV; + (1 — r)C, where z = 1, 2, 3. Since this is an 
affinely symmetric set, we find that the weight (for each point) is one-third the 
area, and r = + 3. Forr = + 3 the evaluation points are the distinct trisection 
points of the median chords and for r = — } the evaluation points are the mid- 
points of the sides. These formulas we consider to offer prospects of extensive 
usefulness, especially the latter. Since the general quadratic function has 6 terms, 
we have used three fewer points than an arbitrary specification of evaluation 
points would have permitted. 
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The cubic polynomial is integrated over the triangle with rV; + (1 — r)C 
and C as evaluation points where r = 3; the weight associated with the centroid 
is (—75)A and the weight with each of the other three points is (25/48) A, where 
A is the area of the triangle. 

The quintic polynomial is integrated precisely with seven points in the triangle 

using rV; + (1 — r)C, weight a; sV; + (1 — s)C, weight b; and C, weight c. We 
hs ee ae a anda = (55 — 18) 4 »-() 
7 , 7 : 1200 ' 1200 : 
c = (9/40)A. Since the general quintic polynomial in two variables has 21 terms 
this formula appears to be a type we call efficient noting that one might not hope 
to accomplish a formula with fewer than 7 = 21/3 points. Here the “3” is the 
number of degrees of freedom for each point due to coordinates and weight. 
However, there are known hyperefficient formulas (cf. [1 ]), which use fewer points 
than indicated by this argument. While we will not reproduce the argument here, 
we used a triangle with vertices (0, 0), (1, —1), (1, 1). Then the requirements of 
the affine symmetry of the formula with the form of the region assured that all 
monomials with odd powers of y could be omitted. This left 12 equations. We 
chose five of these and solved them for a, b, c, r, and s, and verified that the 
remaining 7 were satisfied. 

Over the tetrahedron we have a formula for the general quadratic in three 
variables involving rV; + (1 — r)C, i = 1, 2, 3, 4, and weight a. We calculate 
r = 1/75 and the weight a = (4)A, where A is the volume of the tetrahedron. 
Another formula with points outside the tetrahedron results if r = —1/75. This 
type of formula will generally be less used than one with points inside. 

For the cubic polynomial over the tetrahedron we have r = 3, a = (9/20)A, 
and c = (—#)A; the centroid C must now be included and c is its weight. This 
formula is efficient since it involves 5 points and the number of terms in the general 
cubic is 20. 

Generalization and extension of the affinely symmetric methods is now being 
carried out. However, these specific results seem sufficiently useful to include now. 

4. Conclusion. In this paper we have taken a step towards obtaining reason- 
able numerical integration formulas over simplexes and cones. The cones which 
are not simplexes have not been emphasized. However, the method devised for 
cones permits obtaining a formula for integration over a cone (finite) provided 
a formula is at hand for the base. Thus later on we hope to obtain other formulas 
for the solid sphere which is a special cone based on its surface. 

Error analysis has not been attempted. Experimental calculations on simple 
regions indicated that for “‘reasonable’’ functions there is a decisive factor in favor 
of the classical formulas here presented over Monte Carlo methods. 

P. C. HAMMER 
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Numerical Integration over Simplexes 


1. Introduction. Integration formulae for numerical evaluation of integrals 
over the simplex in m-space have been given inductively by Hammer, Marlowe, 
and Stroud [1] so that it is possible in principle to determine a formula holding 
exactly for the kth degree polynomial in variables. In the same paper certain 
affinely symmetric integration formulae are given for the triangle and tetrahedron. 
Using the theory proposed by Hammer and Wymore [2], it is possible to extend 
the usefulness of methods developed by transformations of the regions and by 
use of Cartesian products. 

In this paper we give two integration formulae of affinely symmetric type for 
the simplex in -space which respectively hold exactly for the quadratic poly- 
nomial and the cubic polynomial in m variables. The method for establishing the 
exact values of integrals needed we believe is new in that the “‘numerical’’ formulae 
are used for the purpose. 

2. The formula for cubic polynomials. Let the vertices of the n-simplex, S,, 


be Vo, ---, V, and then its centroid is given by C = >> Vi/(m + 1). Let A, be 
1 


the hypervolume of S,,. 
THEOREM I: An integration formula exact for the general cubic polynomial over 
S, for n > 1 is given by 











(1) i ft = an © f(Ud) + eof O) 
where 

“ (n + 3)? re —(n + 1)? 

" 4(n + 1)(n + 2) “— 4(n + 2) 
and 

pwn Oho gee oc. 


n+ 3 n+3 


Proof: It may first be remarked that the points U; are on the median lines of 
S, and that the statement of the theorem is in symmetric form. In particular, 
under an affine transformation taking S, onto itself, the set of points {U;} is 
invariant and the centroid C is fixed. Since there exists an affine transformation 
mapping any simplex in E, onto any other we may choose any particular simplex 
S, to carry out the proof. Our choice is specified by vertices as follows: 


Vo= (0, -++,0), VYi= (1, 0, -++,0), 
= (1,1,0, ---,0),*->-, V, = (1,0, ---, 0, 1). 
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It is simply verified that the formula given holds for m = 1 and m = 2. Hence hoi 
we assume that it holds for E,., and proceed to show that it also holds for E, wh 
where m — 1 > 2. Let f be a cubic polynomial in x1, x2, ---,x,. Then f|2,.1 is rg 
a cubic polynomial in x2, ---,x,. Now using a result established in [1] we may out 
write 
. the 
(2) f fdv, = f (f fdt.-1) dx, 
Sn 0 218,-1 : F pol 
= f x1" [an1 Df (x10) + cnsf (x10) Idx: for 
1 n - 
nu: 
where S,_1 is the (m — 1)-simplex with vertices V;, ---, V, in the hyperplane we 
: " 2 n ; thr 
=1,C = i/nd V;, U; = ree ther? ? adhd 1,2, ---, and a,_; and the 
C,—1 are the weights as indicated in the theorem with m replaced by m — 1. It is a 
observed that the hypervolume A,_; is 1/(m — 1)! and A, is 1/n!. Let f be the ’ 
monomial x1‘x2ix3* where 0 < i + 7 + k < 3. Using (2) we find on substitution wit 
and simplification that: wh 
bot 
(3) f seoisestdy Anal (m + 2)?-#-* (37 + 3* + n — 2) — nb] tio: 
x1" ae > * 
* le adie 4(n+i1)mn+itjt+h) Ov 
int 
On the basis of our assumption, (3) gives the value of the integral indicated. On ext 
the other hand, formula (1) applied to f = x,‘x2ix3* gives 
A to. 
4) ————__ { (n+:3)*-*- [nit (n+2)*(31+-3'-—2) J— (n +1) 7}, In 
® aaa mEDy (eta lnt 2)" J-(etaeremy, J In 
Now it may then be directly verified that (4) gives the same result as the right i 
of (3) for O< i +f +k <3, O<i<3 and j >k. Hence in view of the 
symmetry with which the last » — 1 coordinates appear in the set of vertices Uni 
Vo, «++, Va, (4) is verified as correct for all monomials of form x;‘x‘;,x*;, for Mat 
O<itjth <3, i: F ip where i; and iz are taken from 2, ---,. The only 
monomial type thus omitted is x2x3x, provided m > 4. Using formula (2) for this Che 
monomial, we find Alal 
An by t 
5 f x v, = 
6) 4, = Tt ie + 2) + 3) anh 
y 
which coincides with the value obtained on substituting f = xoxsx, in (1). Hence pub 
the formula (1) holds for all required monomials and for the cubic polynomials Mat 
provided it holds on S,_:. Hence by complete induction the formula (1) holds. 
3. A formula for the quadratic polynomial. 
THEOREM 2: The formula 
©) fsa, = “Es 
Sa + n+1°% , 
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holds for the quadratic polynomial over S,, n > 1 provided U; = 1rV; + (1 — r)C 
where either r = 1/Vn + 2 or r = —1/Vn + 2. The choice of the positive sign for 
r gives points U; inside S, for all n. If the negative sign is chosen, the points U; are 
outside S, for n > 2. 

The proof of this theorem may be made along the same lines as that of 
theorem 1, by induction. 

4. Remarks. The formulae given are affinely symmetrical. Those for the cubic 
polynomial are based on rational combinations of the vertices whereas the 
formulae for the quadratic function are based on irrational combinations unless 
n + 2 is a perfect square. The weight of the centroid is negative and increases 
numerically with nm for the cubic case. The formula for the cubic polynomial may 
well be used as an exact integration formula for any polynomial of degree at most 
three. Since the general cubic polynomial has (m + 1)(m + 2)(m + 3)/6 terms, 
the formula (1) is hyperefficient (see [1] or [2]) for > 3 since it is based on 
n + 2 evaluation points. 

The practicality of using simplicial decompositions of regions decreases rapidly 
with increasing m. For n = 1, the formula (1) is based on three evaluation points 
whereas Gauss’ formula uses only two. We conjecture that for no region of 
bounded volume in E, for m > 2 will it be possible to obtain a numerical integra- 
tion formula exact for the cubic polynomial based on fewer than m + 2 points. 
Over the hypercube, Tyler [3] has shown that the cubic polynomial may be 
integrated exactly by a formula using 2” points. Hammer and Wymore have 
extended this result to certain symmetrical regions. 

Extension of affinely symmetric formulae for integration over the simplex 
to higher degree polynomials promises to offer significantly greater complexity. 
In the tetrahedron, for example, we have shown that evaluation points to obtain 
a formula exact for the fourth degree polynomial cannot all be taken on the 
median lines. However, the methods used here may be of use for higher degree 
polynomials. 
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Numerical Quadrature of Fourier Transform 
Integrals 


I. Introduction. In applied mathematical problems the need frequently arises 
to evaluate numerically integrals of the form 


(1a) S(x) = f ” $(k) sin kxdk 
or 
(1b) C(x) = f ” W(k) cos kxdk. 


In some cases the function ¢() or ¥(k) is given by a closed expression which is 
too complicated to permit a sufficiently accurate analytic evaluation of the 
integral for the entire range of the parameter x. In other cases ¢(k) or ¥(k) may 
be available only in numerical form. 

The conventional methods of numerical quadrature (e.g., Simpson’s rule) are 
not suitable for evaluation of the above integrals when x is large. There are two 
reasons for the failure of the standard methods in this case. First, because of the 
rapid oscillation of the trigonometric function when x is large the integrand 
cannot be accurately approximated by simple polynomials unless undesirably 
small intervals of integration are chosen. Secondly, the extremely strong cancella- 
tion between the contributions to the integral from regions where the trigo- 
nometric function is positive and regions where it is negative tends to accentuate 
the errors in the conventional integration procedures. 

Filon [1] and Luke [2] have suggested methods of integration which avoid 
the first difficulty by using polynomial approximations for ¢(%) or ¥(k) rather 
than for the entire integrand. In this note we shall discuss a new scheme which 
attempts to alleviate the second difficulty as well. In this scheme all half cycles 
of the trigonometric functions are treated in an identical manner so that no can- 
cellation errors arise. The integration scheme for the individual half cycles is 
based on a Gaussian type integration procedure [3] which minimizes the error 
in the final result for a given number of integration points per half cycle. In 
developing the integration formulas for the individual half cycles, it is imagined 
that the integrals in (1a) and (1b) are evaluated by first performing the sum of 
the integrand at corresponding points in all half cycles and then integrating the 
result over a single half cycle. The sum over corresponding points in all the half 
cycles possesses certain properties which, particularly for large x, enable the final 
integral to be accurately evaluated with a small number of points per half cycle. 

In the actual numerical evaluation of the integral, the integrations are first 
performed over the individual half cycles, and then the total contributions of 
the various half cycles are summed. Although the integration formulas are not 
accurate for the individual half cycles, most of the error cancels when the sum 
over half cycles is performed (see section V). 

In performing the sum over half cycles, the standard techniques for accelerat- 
ing the convergence of the sums of an oscillating series can be employed. This 
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greatly reduces the number of half cycles which must be considered to obtain a 
final result of specified accuracy. 

A disadvantage of the method is that the values of k at which $(%) or ¥(k) 
must be evaluated depends on the value of the paramter x. Hence the method is 
most useful when results are needed for only a few values of x. If results are needed 
for a large number of values of x, the Gaussian integration feature could be sacri- 
ficed without serious loss in accuracy. This would enable the points at which 
the integrand is to be evaluated to be spaced uniformly instead of at irregular 
intervals as in the Gaussian procedure. It is, of course, still necessary for the in- 
terval in k to be an integral fraction of x/x, but the integral can be evaluated 
for several values of x with the same spacing in k. Only the simplest example of 
the uniform spacing method, i.e., one point per half cycle, will be specifically 
treated in this note. The generalization of the uniform spacing method to a 
greater number of points is, however, straightforward. The discussion in the 
following sections indicates that, although the uniform spacing method is by 
no means as accurate as the Gaussian procedure, the loss in accuracy with a 
corresponding number of integration points may not be serious for practical 
applications. 

It must be recognized that although the method described here is designed 
to reduce cancellation errors associated with the integration formulas themselves, 
the method does not eliminate the possibility of cancellation errors in the nu- 
merical application of the formulas. Therefore, the numerical values of $(%) or 
¥(k) which are used in the evaluation must ordinarily be accurate to a number of 
significant figures equal to the number required in the final answer plus the 
number lost due to cancellation. There are, however, circumstances in which 
consistent approximations may be made in the calculation of ¢(%) or ¥(%) which 
do not introduce as large an error in the final result as would be expected on the 
basis of a consideration of significant figures. These approximations may usually 
be regarded as corresponding to the replacement of the actual physical problem 
by a slightly modified problem which, on physical grounds, must lead to almost 
the same final result. It is, of course, necessary to retain as many significant 
figures in the calculation of ¢(%) or ¥(%) in the modified problem as would be 
necessary on the basis of elementary considerations. This is true despite the fact 
that the integrand as calculated in the modified problem may agree with the 
correct integrand to far fewer significant figures than must be retained in the 
evaluation of the integral. 

In many cases the strong cancellation occurring in integrals of the type appear- 
ing in (1a) and (ib) may be interpreted as resulting from the fact that there is a 
saddle point of the integrand in the complex k-plane which, for large x, is far 
removed from the real axis. Hence the ideal method of numerically evaluating 
the integral would be to integrate along a path in the complex k-plane passing 
through the saddle point. It is interesting to note that one is led to a consideration 
of complex values of k by the following independent argument. One expands the 
function ¢(%) or ¥(&) in a power series in the variable 


(2) u= (1+ RL)", 
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where L? is a suitably chosen parameter. If one then applies the Gaussian inte- 
gration procedure to the infinite integrals in (1a) or (ib) by using polynomials in 
the variable u, one finds that when «x is large the optimum points for evaluation 
of the integrand cluster around the point 7/L in the k-plane. If L has been chosen 
correctly, this point is, in turn, close to the saddle point of the complex integral. 
Although this method of numerically inverting Fourier transforms is in principle 
more powerful than the one considered above, it is probably less practical, mainly 
because of the need to evaluate the integrand for complex values of k. 

Il. Preliminary algebraic manipulation. It is necessary to express the integrals 
(1a) and (1b) as integrations over a single half cycle of the integrand, summed 
then over all half cycles, in order to apply the methods outlined in section I. 
However, in order to simplify the derivation of the method the order of summation 
and integration are interchanged, i.e., the summation is performed first. The 
justification of the interchanging of order of summation and integration is given 
in section V. 

We begin by making the transformation 


(3a) y = kx/x —3 
in (1a), and 
(3b) u = kx/x 


in (ib). For purposes of the derivation, it is convenient to extend the integrals 
to the range — © < k < + o~, which may be done if the definitions of the func- 
tions @ and y are extended to negative values of k by the definitions 


(4a) ¢(—k) = — $(), 

and 

(4b) v(-k) =v). 

In most applications of interest, the functions thus defined are regular in the 
range k = — © to +. (Note that the restrictions implied on ¥(k) and ¢(k), 


namely ¢(0) = 0 and y’(0) = 0 are physically reasonable. For example, if (1a) 
arose from the inversion of a three-dimensional Fourier transform f(k), then 
o(k) ~ kf(k).) 

Applying (3) and (4) to (1) leads to 


Tr 4 
(Sa) Ste) = = J cos ayo, 2)dy, 
and 
(5b) C(x) = ~ f cos wny (x, x)du, 
where 
+00 7 
(6a) co.2) = E (-)(Z+n43)), 


n=—o 
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and 


+o . 
(6b) vu,2) = © (-¥(Ztu +2). 


The function o(y, x) has the following properties [all the subsequent remarks 
apply equally well to y(y, x)] as may be seen from (4) (6): 


(a) o(y, x) = o(—y,x), 

(b) o(3, x) bas o(—}, x) =0 

(c) oly +n, x) _ (—)*e(y, x), 

(d) ¢ is regular, —$ < y < $ if $(&) is regular for -~> <k< +a. 


Because of these properties, ¢ (and y) can be expanded in Fourier series of 
the form 


(7a) o(y,x) = > a,(x) cos (2m + 1)xy, 
n=0 
or, since cos ”@ is a polynomial of exact degree in cos @, 


(7b) o(y, x) = cos ry > a, (x) cos ay, 


with a similar expression for y. 
Ill. Gaussian integration procedure. Substitution of (7b) in (5a) gives 


3 wo 
(8) S(x) = — cos’ ry >> an(x) cos wrydy. 
2x 4 n=0 


In order to integrate such an integral as accurately as possible, the Gaussian 
integration procedure is employed [3]. A set of polynomials in cos* xy orthogonal 
on the interval (—}, 4) with weighting function cos’ ry is found. Representing 
these polynomials by I,,(cos zy), one has 


4 
(9) f cos? ryI', (cos ry)I',(cos ry)dy = 0 for nm, 
ae | 


The polynomials in cos? xy which obey (9) may easily be shown to be related to 
the Chebyshev polynomials of the first kind [4], 7,(x): 


(10) T,.(cos ry) = (cos ry) T2n41(cos ry) 


= (cos ry)! cos (2m + 1)ry. 


See, for example, [4]. (The authors are indebted to the referee for pointing out 
this relationship.) The 2N-point Gaussian quadrature formula is simply 


ry N 2wW™ om 
(11) [ics ryo(y, x)dy = dos ays 70% »*)s 
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where the y;™ are the zeros of I'y(cos ry) : 


4%j-1 
LID apn aint = ane N. 
(12) 9 2(2N 4 1) ’ J 1, 2, ’ ’ 


and the W, are the so-called Christoffel numbers. (Since o(y, x) is an even function 
in y, the contribution to the numerical integration from —+; is equal to the con- 
tribution from +~;; furthermore, it is easy to show that the Christoffel number 
associated with +-y; is equal to the one associated with —+,;. For this reason, we 
consider only the positive zeros of 'y(cos zy), and introduce the factor of two on 
the right side of (11). Thus, although the integrand is evaluated at only N points, 
the formula is called a 2N-point formula. In section V, however, where the order 
of integration and summation is reversed for numerical convenience, it is neces- 
sary to evaluate the integrand at both ++;, which gives the numerical integration 
formula as stated there a slightly different appearance.) 

Note the factor of (cos ry;™)- in (11), which arises from the fact that a 
cos ry was factored out of o(y, x) in (8), and the weighting factor cos* ry was 
used, although the original integrand involved only the first power of cos zy. 

The W;™ may be found by requiring that (11) give an exact answer if 


a(y, x) is a polynomial of degree up to and including N — 1 in cos* ry. It then 
cos ry 


follows automatically that (11) will also give an exact result for polynomials 


from degree N to 2N — 1 inclusive in cos? ry. Since 


“2 1 TC 
as [loot 99 = era 


then the W;™ are the solutions of the NV simultaneous equations: 


Be: Ls eae pee Lak o 
(14) Vrad¢i 22% (fe + ) ah lite 
A= 1,2, +++, MN 


(It is also possible to obtain a general closed expression for the Christoffel num- 
bers [3], although the actual evaluations of the Christoffel numbers are just as 
simple if (14) is used directly.) 

The y; and W; are given below for the two- and four-point formulas: 


(a) Two-point Formula (N = 1): 


nm=% 
cos ry; = $3 
Wi = 3 
(b) Four-point Formula (NV = 2): 
y= Yo, y2 = ts 
cos ry: = 0.58778 52 cos ry2 = 0.95105 65 


W, = 0.06909 832, W2 = 0.18090 169. 
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It will sometimes be found that a one-point formula gives sufficient accuracy 


(see section IV). This formula gives an exact answer only if eG 





= is independent 
of y. In this special case, (11) becomes 


3 
(15) £ cos rya(y, x)dy = $e(0, x). 


IV. Choice of N. The number of points which must be used in the Gaussian 
procedure in order to obtain a numerical answer within a required degreefof 
accuracy may be shown to depend upon the magnitude of x. From the arguments 
given here, a qualitative dependence of the error introduced by the}N-point 
formula may be obtained as a function of x. Consider the example 


(16) ¥(k) = (1 + LR). 


The integral (1b) may be evaluated analytically, yielding 
‘-_. in 
(17) C(x) = aL ee, x x/L. 


The 2N-point quadrature formula gives for C(x), 





(18) Cw) = FE O™.2) 
where 

+00 x2 
(19a) v(yj,x)= LD (-)* 





x? + L¥x*(yj™ + nn)? 
x* sinh x* cos ry;™ 
sin? ry; cosh? x* + cos? ry; sinh? x*’ 





1 (95, x) », x) 


(19b) cose = 2x*[e* — e-*™*][1 + (4 cos* xy — 2) 4 --- J, 


Thus, it is seen from (19b) that for large x, (7) is rapidly convergent. From (18) 
and (19b), one sees that in the limit of x large 


(20) C(x) == oe + O[e-***], 


independent of N, so that even the one-point formula may be used. As x decreases, 
however, it is necessary to go to larger values of N. 

In Figure 1, the error introduced into the evaluation of C(x) by the one-, two-, 
and three-point formulas with ¥(%) given by (16) is plotted as a function of x*, 
where the numerical evaluation of C(x) is obtained from (18) in conjunction with 
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(17a). Note that the error depends on the ratio x/L (i.e., x*) since this quantity 
determines how much the integrand decreases in one cycle of the trigonometric 
factor. When the rate of decrease is small (i.e., x* is large) the required number 
of integration points per half cycle is small. 

V. Outline of numerical procedure and discussion. In order to evaluate the 
sum appearing in (11), it is convenient to change the order of summation and 
{numerical) integration. That is, the order of the two sums on the right side of 


10> 
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Fic. 1. Per Cent Error in Integration Schemes as a Function of x*. 





(11) are reversed. This procedure corresponds to a separate integration of every 
half cycle, with a subsequent summation over the cycles, as is described in section 
I. It may easily be shown that this reversal may be made since if the integrals 
(1) exist, then the series (6) for o(y) and y(y) converge uniformly on the range 
—} <y < }. (The only restriction on ¢(k) (and ¥(k)) is that there exist a ko 
such that $(k:) > ¢(k2) if ko < ki < ke.) It is also convenient to change the 
formulation slightly so that the half-cycle sums are carried out from zero to 
infinity. 
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ntity We define 
etric +s 
mber (21a) S(x) ==¥ S,(x), 
n=0 
a diies and 
d « 
le of mn C(x) = *[ scat +z cute) | 
with 


(2a) S42) = (- fo] 20 +0 + ¥| cos ayy 


= (-po 


jai COS wy 


x {o(Zo +0443) +6(Zt- 9 +0449), 


and 


(22b) Ca(e) = (—)" fv (Zu + 2]) 008 xudu 


=(-y 5 


x fy (E tum + 0) +¥(EC- a +23)}- 
t x x 


In order to sum the series in (21) any standard method [5 ] [6] for the summa- 
tion of oscillating series may be employed. A method suggested to the authors 
by Mrs. M. Ray has been found to be quite convenient; it appears to be suffi- 
ciently accurate for most purposes, and is quite simple. The method is to form 
the sequence of partial sums, 


(23) B, = > Sa(x), 
a=0 

so that 

(24) S(x) = = Lim Bu. 

The limit of the sequence fo, 61, ---, 8, °-* is then found by the “‘averaging”’ 
ed technique. That is, we define 
ion 
rals (25) ai = (Bn + Bn41), 
nge a, = (a, a a"), 
i Ro etc., 
the , ; 
- to (26) Lim 6, = Lim qj. 


no reo 
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The sequence of a;’s in general converges much more rapidly than does the 
sequence of 8,’s. (Note that a;, as defined by (25), is given in [7] (Exercise 120, 
p. 271), as the j-th partial sum of a series equivalent to the original oscillating 
series. Thus, our procedure corresponds to summing the equivalent series which 
is defined in [7] (paragraph 144, p. 244). The summation procedure described 
here is useful because the equivalent series converges more rapidly than the 
original series. Still faster convergence could presumably be obtained in typical 
cases by the application of the more refined procedures discussed in the references 
given in [6] and [7].) 


TABLE I. Replica of Portion of Work Sheet, Demonstrating a Convenient 
Method of Setting Up Calculations 


m o(Ztit+n+4]) S.=(-»We(ZDit+e+4]) B= 25, 


0.27763 71 


+0.13881 86 


+0.06689 27 


a; aed 4(Bn + Bn+1) 


+2.6191 X 10-3 
il 0.25709 43 —0.12854 72 —0.06165 45 —1.8693 x 10° 
12 0.23914 07 +0.11957 04 +0.05791 59 +2.0741 x 10-3 
13 0.22336 71 —0.11168 36 —0.05376 77 —1.4098 x 10° 
14 0.20943 14 +0.10471 57 +0.05094 80 +1.6853 X 10-% 
15 0.19705 09 —0.09852 545 —0.04757 74  —1.0791 X 10° 
16 0.18599 31 +0.09299 66 +0.04541 92 





Table I is the replica of a portion of a work sheet in which the integral 
(27) S(x) = Cx i+Fe —— sin kxdk = }xe7*, 


was evaluated by the single-point formula, with x = 10, demonstrating a con- 
venient method of setting up the calculation. In Table II, the summation process 
is illustrated for this case. The criterion for ending the calculation is the near con- 
stancy of the last two or more terms in the upper diagonal of Table II (a{” = a$"),) 


TABLE II. Illustration of Summation Procedure 


n a™ X10 a” X10 a3 K 108 ay K 104 as™ X 104 
10 
+2.6191 
11 3.749 
— 1.8693 2.386 
12 1.024 2.271 
+2.0741 2.173 2.270 
13 3.322 2.270 
— 1.4098 2.350 
14 1.378 
+1.6853 2.204 
15 3.031 
— 1.0791 
16 


which indicates the sequence has closely approached its limit. We then take 


S(x) =- = ay. 
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The numerical result for S(x) is 7.131 X 10~-*, compared to the rigorous value 
of 7.1314 X 10-5. Since the error introduced by the quadrature rule is only of the 
order e~** = ¢~**, the main error in this procedure is introduced in the summation 
process, during which significant figures are lost. Thus, the numerical result ob- 
tained here agrees with the analytical result to as many significant figures as were 
retained during the numerical procedure. However, three of the seven figures to 
which ¢(k) was originally evaluated were lost during integration. 

The same integral has been treated numerically by Simpson's and Filon’s 
methods, and it has been found that in order to obtain accuracy comparable to 
that obtainable by the method described here the integrand must be evaluated 
at at least 10 times as many points. 

The tables and graph in this paper were prepared by Miss D. M. Keaveney. 

H. Hurwitz, Jr. 


P. F. ZWEIFEL 
General Electric Co. 
Knolls Atomic Power Laboratory 
Schenectady, New York 
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Formulas for the Partial Summation of Series 


1. Introduction. The paper gives a table of the coefficients A,,(m) in the 
10 

“partial” summation formula S, = >} Aw»(m)Sn. The coefficients A,,(n), 
m=4 


m = 4(1)10, are tabulated exactly in the fractional form C,,(m)/D(n) for 
nm = 11(1)50(5)100(10)200(50)500(100)1000. For every n, except 47, D(n) is the 
least integer containing no more than ten digits exclusive of final zeros, to permit 


10 
ready division on a ten-bank calculating machine using S, = [ > Ca(m S| / D(n). 
m=4 


The purpose of A,,(m) is the calculation of the sum of m terms of a slowly con- 
vergent series, or, more generally, the evaluation of the m-th term of a sequence S,, 
which is either slowly convergent or asymptotically characterized by S,, ~ f(m) 
(convergent or divergent), in such a manner that the auxiliary sequence S,,/f(m) 
is slowly convergent. The formula for S, was obtained by Lagrangian extrapola- 
tion upon 5S, considered as a polynomial in 1/m, based upon the last 7 values of 
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the first 10 terms of the sequence S,,. Examples of the calculation of S, from the 
initial values S,, Ss, ---, S10 are given for S, = > 1/r?, } 1/r, the smallest 
r=1 r=i 


zero of the m-th order Laguerre polynomial, and the m-th zero of the Bessel 
function Jo(x), all the results being of high accuracy. 

2. Relation to earlier work. In a previous article [1] the writer gave a 
method, with formulas, for the summation of a slowly convergent series or, more 
generally speaking, for the evaluation of the limit of a slowly convergent sequence. 
Those formulas were concerned only with the ‘complete’ sum or limit, i.e., 
limit S,,. But in many practical problems we wish to have S,, for some larger finite 


mo 

value of m, say m = n. For example we may need to estimate the sum of 100 or 
more terms of a series from the sums of the first few. Then we may have to know 
the smallest zero of a sequence of polynomials or transcendental functions for 
some large order from a tabulation of those zeros for only lower orders. The pro- 
duction of this newer table of formulas for “‘partial’’ or “‘incomplete’’ summation 
was motivated by the extreme accuracy of the method in the preceding tables of 
formulas for complete summation, namely, Lagrangian extrapolation for S,, con- 


1 1 
sidered as a polynomial in ny Only here, instead of extrapolation for < 0 


(m = «) from the last tabulated values up to a fixed m, we extrapolate for 


| ee | : : ‘ } 
— = — (m = n) when we require the value of S,. But in employing this method, 
m n 


the sequence S,, must always be examined, usually in the roughest way, to see 
whether it converges slowly as it is, or whether it diverges in some regular manner, 
in which case it is essential to calculate with some auxiliary sequence obtained 
by dividing through by some simple function of m, often m itself, to insure con- 
vergence. Even when the sequence S,, does not diverge, knowledge of its asymp- 
totic behavior, say S, ~ f(m) as m— ©, enables one to obtain much more 
accurate results by extrapolating upon the auxiliary sequence S,,/f(m) instead 
of S,, (as will be apparent from some of the examples below). 

3. Choice of arguments. The earlier article [1], which was devoted only to 
complete summation, gave sets of multipliers depending upon the number of 
terms of S,, that were available, and also the degree of accuracy that was desired. 
But since ‘‘incomplete”’ summation requires a separate set of multipliers for every 
desired S,, where m should assume a number of values to be generally useful, to 
save excessive space it was decided to choose a single number of available values 
of S,, as well as just one degree of approximation, which would serve best in an 
all-purpose capacity. Careful scrutiny of the examples for complete summation 
in the previous article [1] showed that the best single all-purpose choice of the 
total number of available values of S,, and the number of end values upon which 
to base the approximate formulas was 10 and 7, respectively. Other considerations 
for the choice of 10 and 7, besides adequacy of the approximation, are: 

A. One does not often go beyond the computation of the first 10 terms of a 
sequence (especially when the calculation may be involved, as in the computation 
of the roots of transcendental equations) and 

B. From the examples in [1] where the results of the 7-point formula were 
comparable with those of the 4-point formula, the substantially greater accuracy 
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in the former justifies the somewhat larger table of formulas that is needed for 
7 instead of 4 points. 

4. Range of table. This present table gives the coefficients A,,(m), more 
simply written as A,,, the m being understood, in the following summation 
formula: 


(1) S. =F An(n)Sa. 


The coefficients A,,, m = 4(1)10, are tabulated exactly in the fractional form 
Cn(n)/D(n), more simply written as C,,/D, the m being understood, for m = 11(1) 
50(5)100(10)200(50)500(100)1000. The left-to-right direction of the columns of 
Cn as m goes from 10 down to 4 was chosen to comply with the natural tendency 
in extrapolation to start the multiplication with the last or 10th term of Sa and 
to proceed backwards. The denominator D(m) = D was chosen as the least 
common denominator of the A,,(m), and nearly every D had no more than ten 
digits exclusive of final zeros, so that use of (1) in a ten-bank desk calculating 
machine in the convenient form 


(1’) S, = (1/D) r CouSim 


permits ready division by D. In the few cases where the l.c.d. of the A,,(m) had 
more than ten digits for division, both numerators and denominator were multi- 
plied by 2, 4, or 8, resulting in larger values of C,, and D, but with no more than 
ten digits in D exclusive of final zeros, with the sole exception of m = 47, where 
division by D = 47® is conveniently managed by two successive divisions by 
47% = 103823. 

At the suggestion of the referee, decimal values of A, have been computed 
and listed by R. B. Horgan [4]. 

5. Error terms. In addition to the theoretical error in formula (1) or (1’), the 
user should always bear in mind the size of the coefficients A,,(m) and the fact 
that in most practical examples the quantities S,, are inexact, so that if the error 
om in S» satisfies | om| <e, an upper bound for the error in (1) or (1’) besides the 

1 
theoretical one, is e > |A,|. In many cases the actual error due to in exact 
4 


values of S,, is a small fraction of that upper bound. 

6. General formula for A,,(n). The arguments 2 were chosen to fulfill most 
of the needs in summation or evaluation of sequences. For all other values of 2 
the user may calculate A,,(”) from the formula 


10 
II ( — j) 
m*® pel 
2 An oe Dd : ? ~ 0 
” = 58 G@— mil m— + 
where j = m is absent from J]’. 
7. Method of computation. The computation of this table used (2), the pro- 
cedure being to first get the seven multipliers m*/T]’(m — j),m = 4(1)10, which are 
10 





independent of m, and to multiply each by the corresponding [] (”—j)/(n—m)n*. 
j=4 
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For the convenience of the computer, we list those multipliers: 











m 10 9 8 7 6 5 4 
m’ 12500 | —17714 7 | 16384 | —11764 9 | 972 | —3125 | 256 
Wm-j| 9 40 3 36 i 24 ~«| «45 


A preliminary functional check upon the A,,(™) was performed by using 





10 

(3) >> An(n) me 1, 
m=4 

and a final functional check made use of the formula 
10 An(n) 1 

(4) m> _—s n 


The author is grateful to Miss Isabelle Arsham for her assistance in the calcu- 
lation and checking of the table of A,, (7). 

8. Illustrations. The following examples are intended to illustrate the wide 
applicability and great accuracy of this table of A,,(m) for many different types 
of problems: 

Example 1: As the most elementary type of problem consider the sum 
Sn = > 1/r*, where the 4th to 10th partial sums are, to 10D: 


r=1 


S, = 1.42361 11111 
Ss = 1.46361 11111 
Ss = 1.49138 88889 
S7 = 1.51179 70522 . 
Ss = 1.52742 20522 
So = 1.53976 77312 
Sio = 1.54976 77312 


Suppose that we want both S;5 and Sx, and employ (1), using the table of A,,(”) 
form = 15 and m = 24. The calculated value of S;; is 1.58044 02838, which agrees 
with the true value of Sis = 1.58044 02834 to around 4 units in the 10th decimal. 
The calculated value of Sx is 1.60412 34045, which agrees with the true value 
of Se = 1.60412 34036 to within a unit in the 9th decimal. 


Example 2: Consider the partial sums of the divergent series S, = >> 1/r, 


r=1 


whose values from m = 4 up to m = 10 are, to 10D: 


Ss 2.08333 33333 
Ss 2.28333 33333 
Se = 2.45000 00000 
S; = 2.59285 71429 
Ss = 2.71785 71429 
So = 2.82896 82540 
Sio = 2.92896 82540 
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Suppose that we wish to calculate S2., whose true value is 3.85441 97162. Al- 
though S,, diverges, the sequence s,, = S, — log m converges slowly to Euler’s 
constant, so that one bases the use of (1) upon s,,, m = 4(1)10, to calculate so. 
and to obtain Sy from sex + log 26: 


$4 = 0.69703 89722 
Ss = 0.67389 54209 
Se = 0.65824 05308 
$7 = 0.64694 69939 
Ss = 0.63841 56012 
Sp = 0.63174 36767 
Sio = 0.62638 31610 


Employing s, in (1), we find for = 26, sxx = 0.59632 31116. Since log 26 
= 3.25809 65380, Ses = ses + log 26 is found to be 3.85441 96496, which agrees 
with the true value to 3 of a unit in the 8th significant figure. 

Example 3: From the known values of the smallest zero of the first ten 
Laguerre polynomials L,,(x), one wishes to calculate the smallest zero of the 
fifteenth Laguerre polynomial [2]. It is true that the »-th zero of L,,(x), or x,™, 
converges to zero as m— © in accordance with the asymptotic approximation 
2(x,™)* = m—+(vr + 0(1)). But a much smoother approximation can be had from 
the sequence S,, = mx,™; these values from m = 4 to 10 are, to 9D: 


S, = 1.2901 90758 
Ss = 1.3178 01599 
Ss = 1.3370 79625 
S; = 1.3513 05736 
Ss = 1.3622 37058 
Ss = 1.3709 00049 
Sio = 1.3779 34705 


Use of (1) for m = 15 yields Sis = 1.3996 17169, from which S,5/15 = 0.09330 
78113, agreeing with the true value of x,;“ = 0.09330 78120 to within a unit in 
the 9th decimal place. 

Example 4: From the first ten tabulated zeros of the Bessel function Jo(x) [3], 
to calculate both the 11th zero and the 42nd zero. The sequence Jo,.. of zeros 
of Jo(x) diverges, but (from the asymptotic formula), jo../m converges, and 
Sm = jo,m/m yields a good approximation ; the values for m = 4 to m = 10, are, 
to 9D: 

S, = 2.9478 83610 
Ss = 2.9861 83542 
Ss = 3.0118 43995 
Sz = 3.0302 33804 
Ss = 3.0440 58941 
Sp = 3.0548 31015 
Sio = 3.0634 60647 
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22 
Bas 
HEE 


2 10141 25000 
162 31600 00000 
60868 50000 


232 47840 00000 
20540 78125 
326 26230 00000 
23 98987 50000 
29 97592 00000 
32 78616 25000 
609 64540 00000 
44 11907 50000 
116 35800 00000 
11 70852 37500 


72 03784 00000 
90190 10000 

2221 01488 00000 
4 05619 50000 
5728 75996 00000 
54 65473 50000 
2596 32172 80000 
10 88321 93750 
53325 11184 00000 
289 70259 50000 


426 47618 65000 
1854 50345 50000 
246 25783 87500 
1403 31064 50000 
2525 53336 75000 
4540 42488 90000 
33154 45489 12500 
1578 07226 50000 
66264 52127 75000 
11388 63953 50000 


7 23664 40627 25000 
11 09639 82405 50000 


8314 03097 25000 


64 44977 33995 50000 
352 28055 64733 00000 
125 38767 78511 50000 


42 15203 66228 50000 
275 89915 32841 50000 
4333 64236 10980 50000 
20 44825 96356 50000 
4173 69687 28899 50000 


—Co(n) 
111 
37 
2008 
47 
818 
1578 
47877 
1773 
1 43632 
1450 


51677 
17075 
03347 
12015 
24428 
69974 
10 13339 

36302 
53 64737 

17479 


90 09678 
7 18162 
16 12843 
1 
1 


w 


36002 
28996 
48307 
340 55121 

6 46960 
500 31598 
1 87827 


718 12596 
63512 
1009 73689 
74 30964 
92 92752 
101 71796 
1892 78696 
137 07292 
361 74885 
36 42374 


224 70962 

2 81945 

6955 56538 
12 72174 
17990 21499 
171 82003 
8169 81995 
34 27430 

1 68057 75451 
913 60782 


2 10126 72999 
36124 21246 


22 97934 88571 
35 26089 65656 
26432 75819 

204 98191 74911 
1120 75082 63074 
399 00262 63551 


134 18049 83957 
878 47158 42992 
13800 98178 29522 
65 12908 25497 
13295 01253 33637 


60261 


12417 
46949 


27124 
34607 
78739 
47797 
81512 
63831 


52039 
69631 
42307 
38921 
91153 


Cs(n) 


91 75040 

34 40640 

1981 80864 
49 15200 

865 07520 
1703 11680 
52481 22880 
1968 04608 

1 61021 95200 
1640 03840 


58825 11360 
19552 66560 
9 24623 83104 
13891 58400 
3 76543 64160 
81497 29280 
11 83855 41120 
42529 25952 
63 00631 04000 
20575 02720 


106 27448 83200 
8 48719 87200 
19 09362 52416 
13 47010 56000 
1 53183 84640 
57443 94240 
405 47804 77440 
7 71221 74976 
597 07490 30400 
2 24385 63840 


858 73970 38080 
76018 48320 

1209 60610 46784 
89 09135 87200 
111 49875 60960 
122 13515 05920 
2274 28873 01120 
164 80909 88544 
435 21933 31200 
43 84741 78560 


271 21051 23840 
3 41001 83040 
8426 95127 85920 
15 43510 42560 
21853 97089 07520 
208 94140 82560 
9943 97965 51680 
41 75079 01440 


2 04862 57586 99520 


sn = 


28 
43 


252 
1378 
491 


165 
1081 


1114 39267 43040 


1644 09078 57920 
7161 98637 15840 

952 45103 92320 
5434 43528 90880 
9790 93084 56960 
17618 76032 71680 
28759 13219 27680 
6133 07876 96640 
57699 37361 30560 
44315 47415 85920 


22047 17200 17920 
33395 84199 06560 
32500 91815 73120 
13382 69862 29760 
95553 32454 80960 
04016 57646 28480 


18859 10014 77120 
74412 43622 60480 


16997 54339 11945 62560 


22554 64819 91680 


16378 55590 39777 38240 





TABLE OF 


—Cr(n) 
41 
16 
988 
25 
452 
905 
28263 
1070 
88324 
905 


32686 
10920 
18708 
7823 
12803 
46200 
72991 
24237 
35 98882 
11776 


60 94424 
4 87553 
10 98606 
7 76189 
88390 
33188 
234 54768 
4 46607 
346 12071 
1 30202 


498 75081 
44189 

703 71461 
51 87085 
64 96436 
71 21093 
1326 89717 
96 21535 
254 23175 
25 62783 


158 90962 

2 00203 

4955 69306 

9 08966 

12884 96683 
123 31650 
5874 12798 
24 68240 

1 21195 43569 
659 67523 


974 26207 
4247 76149 
565 30709 
3227 48153 
5817 86105 
10474 04799 
76575 96983 
3648 78220 

1 53363 04256 
26380 70362 


16 81746 04905 

25 84226 93115 

19391 64584 

150 49140 16471 
823 29572 38353 
293 23935 52389 


98 68113 06717 
646 37486 41980 
10158 39760 93632 
47 95270 76172 
9790 95570 52294 


wn 


nN 


a 


57951 
22330 
77800 
46170 


19655 
30330 
17040 
25820 
37390 
52770 
03360 
72530 
16940 
74030 


59590 
86730 
21190 
82810 
19470 
54110 
55960 
39130 
60970 
73690 


78440 
84670 
02590 
56570 
04720 
40110 


90590 
26110 
11670 
09010 
74930 








\BLE OF 


) 


41 17715 
16 47086 
88 25160 
25 21050 
[52 94865 
0S 89730 
163 99576 
170 60590 
$24 98675 
0S 89730 


586 25360 
20 18018 
108 SS855 
323 65850 
303 51120 
200 76230 
91 10417 
137 37470 
382 91000 
176 66490 


424 08575 
553 92686 
506 36200 
189 27750 
390 53405 
188 78290 
168 17376 
507 36890 
971 08975 
202 14830 


81 74880 
189 30054 
461 10995 
B85 58550 
436 60120 
093 96670 
717 57951 
535 22330 
175 77800 
783 46170 


962 19655 
203 30330 
306 17040 
966 25820 
683 37390 
650 52770 
798 03360 
240 72530 
569 16940 
523 74030 


207 59590 
149 86730 
709 21190 
153 82810 
105 19470 
799 54110 
983 55960 
220 39130 
256 60970 
362 73690 


905 78440 
115 84670 
584 02590 
471 56570 
353 04720 
389 40110 


717 90590 
980 26110 
632 11670 
172 09010 
294 74930 
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Am(n) = Cn(n)/D(n) 
Ce(n) 


9 79776 

4 08240 

251 94240 
6 56100 
119 75040 
242 49456 
7642 25280 
291 89160 
24249 45600 
250 19280 


9073 65888 


99 82611 07200 
37586 65680 


144 10366 24896 
12778 00920 
64693 


15 02180 31600 

18 82672 24320 

20 65056 11676 

385 02943 31520 
27 93569 40720 

73 85691 45600 

7 44917 56920 


46 29945 48480 
58443 18480 

1448 97504 99840 
2 66126 95395 
3776 79708 80640 
36 18202 54320 
1725 00930 20160 
7 25382 11880 
35641 78576 35840 
194 11799 75280 


286 98566 76540 
1252 31690 15280 
166 78169 60040 
952 77671 51760 
1718 37290 32470 
3095 03664 38160 
22636 89142 54560 
1079 00963 20080 
45366 41990 62620 
7805 89094 59440 


4 98146 85057 61440 

7 66004 75168 47920 
5750 83455 86340 

44 64661 54117 18320 
244 31894 33391 90720 
87 04077 55242 93360 


29 30103 02575 41840 
191 97238 00716 45360 
3017 57739 70941 79920 
14 24648 08957 91760 
2909 16581 42133 49680 


—Cs(n) 


1 09375 
46875 

29 $3125 
7812S 

14 43750 
29 53125 
938 43750 
36 09375 
3016 40625 
31 28125 


1139 53125 
383 90625 
18370 62500 
278 90625 
7630 87500 
1665 46875 
24376 40625 


= 
~ 
= 
o 
x 
a 


88 
31441 40625 
431 68125 


2 24142 18750 
17986 71875 


_ 


40645 
28792 96875 
3286 96875 
1237 03125 
76103 59375 
16715 78125 
12 97931 25000 
4891 21875 


18 76778 75000 
1665 


56247 34375 
96071 09375 
45896 21875 


69886 

35017 81250 
65521 40625 
66895 31250 
97571 03125 


6 07819 27500 
7686 65625 

190 86846 62500 
35101 68750 
498 70901 43750 
78228 93125 
228 19233 93750 
96028 40625 

4721 49427 75000 
25 73016 46875 


38 07823 09375 
166 30058 15625 
22 16320 54875 
126 68776 65625 
228 60431 34375 
411 93370 96875 
3014 04135 37500 
143 71729 55625 
6044 39890 03125 
1040 30841 90625 


66458 97608 62500 

1 02265 53463 21875 

768 14416 59375 

5 96568 31390 40625 

32 65522 01225 25000 
11 63635 64672 46875 


3 91854 54213 46875 
25 67946 25583 71875 
403 72375 76966 34375 
1 90631 47076 53125 
389 31783 01659 28125 


ewSnne® 
3 
w 
x 
wn 


— 


Ca(n) 


37117 62432 


é 
: 


83808 

20246 28224 
1510 37952 
48204 29824 
20619 21792 
05536 


4146 75968 
01836 19584 
11245 79648 


64796 83584 
20314 72128 
96064 01024 
49437 39136 
91946 47552 
88228 87936 
89219 01056 
24342 77248 
66401 20832 
45 21987 15904 


2891 83356 88704 
4452 95462 36672 
33 46373 43744 
25998 64798 46912 

1 42352 93948 39552 
50737 51712 78976 


17091 66070 20544 
1 12034 08376 11776 
17 61680 06252 95872 
8319 52049 17616 
16 99247 76781 56288 


» = 3 os 
- 
SasSicu -_— oso = @ 


D(n) 


17 71561 
1 86624 

48 26809 
67228 

7 59375 

10 48576 
241 37569 
7 08588 

470 45881 
» 4 00000 


122 52303 
35 43122 
1480 35889 
19 90656 
488 28125 
96 53618 
1291 40163 
43 02592 
5948 23321 
18 22500 


8875 03681 
671 08864 


35187 43761 
128 00000 


47501 04241 
40 84101 
63213 63049 


19773 26743 
1953 12500 


11072 25625 
129 60000 

3 01675 56250 
$25 21875 

7 11914 06250 
6553 60000 

3 01719 61250 
1230 18750 

58 80735 12500 
31250 00000 


44289 02500 

1 86624 00000 

24134 04500 

1 34456 00000 

2 37304 68750 

4 19430 40000 

30 17196 12500 
1 41717 60000 

58 80735 12500 
10 00000 00000 


610 35156 25000 


27679 21875 00000 
9765 62500 00000 


3240 00000 00000 
21008 75000 00000 
3 27680 00000 00000 
1537 73437 50000 
3 12500 00000 00000 
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Use of the above values of S,, in (1) for m = 42 gave S42 = 3.1229 14703, from 
which 42 X Ss = 131.16241 75, agreeing with the true value of 131.16244 63 
to within 3 units in the 8th significant figure. But the same formula (1) for the 
very next or 11th zero was so accurate as to give an answer of 33.77582 019, which 
agrees with the true value of 33.77582 021 to 2 units in the 10th significant figure. 


HERBERT E. SALZER 
Convair 
San Diego 12, California 


The work reported was done while the author worked at the Diamond Ordnance Fuze Labora- 
tories, Washington, D. C. 
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Table of Coefficients for the Partial Summation of Series 


This note supplies a table of the coefficients A,,(m) as described by Herbert 
E. Salzer [1]. The computation of the table was undertaken at the suggestion of 
the referee of Dr. Salzer’s paper and the Chairman of the Editorial Committee. 
Responsibility for its correctness does not fall upon Dr. Salzer. 

For m = 4(1)10 and m = 11(1)50(5)100(10)200(50)500(100)1000, the table 
lists 


m* TT (nw — 3) 
”) An(™) = Gu — m) I Gm —D) 





to fifteen decimals. (j = m is absent from J]’.) 
The computation was performed on the SWAC machine by means of double 
precision routines. The seven values of 


més 


a Tos 


were stored as constants in the computing routine. For each m, the seven values 


II (@ — 4) 
of omg S were obtained and then multiplied by (2), resulting in (1). 
n’(n — 


The results were converted to 25-decimal numbers and punched on cards by 
SWAC. They were listed on an International Business Machines Corporation 











TECHNICAL NOTES AND SHORT PAPERS 


Coefficients for Partial Summation of Series 


A,,(n) 

3.95131 75103 76442 

9.37714 33470 50754 
17.40280 17267 72284 
27.89016 48122 80597 
40.55967 07818 93004 
55.07469 17724 60937 
71.09249 48572 90725 
88.29178 59179 09984 
106.38550 90735 78621 
125.12500 00000 00000 


144.29940 23246 07871 
163.73271 93362 23815 
183.27988 01917 55393 
202.82258 71270 57613 
222.26534 40000 00000 
241.53185 36532 10640 
260.56185 16990 72116 
279.30837 96929 85066 
297.73546 82433 50869 
315.81618 65569 27298 


333.53101 10110 96866 
350.86646 67606 35376 
367.81400 03486 21887 
384.36904 62365 94911 
400.53025 52507 88361 
416.29885 77283 27321 
431.67813 84464 44154 
446.67300 41679 95069 
461.28962 78468 17233 
475.53515 62500 00000 


489.41747 00281 06694 
502.94498 71587 40688 
516.12650 22606 36098 
528.97105 55761 27495 
541.48782 64943 80373 
553.68604 73746 33594 
565.57493 41604 04691 
577.16363 08793 01430 
588.46116 56213 25944 
599.47641 60000 00000 


650.61572 25181 63361 
695.91126 54320 98766 
736.22631 59781 13077 
772.28678 52680 43077 
804.69824 40384 08780 
833.96507 26318 35937 
860.50810 76723 17788 
884.67972 36193 66967 
906.77629 08297 96555 
927.04830 40000 00000 





m 
10 
10 
10 
10 
10 
10 
10 
10 
10 
10 


10 
10 
10 
10 
10 
10 
10 
10 
10 
10 


— 
Oo 
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n 

110 
120 
130 
140 
150 
160 
170 
180 
190 
200 


250 
300 
350 
400 
450 
500 
600 
700 
800 
900 


1000 


11 
12 
13 
14 
15 
16 
17 
18 
19 
20 


21 
22 
23 
24 
25 
26 
27 
28 
29 
30 


31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
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A,,(n) 


962.93875 

993.71112 
1020.37531 
1043.69507 
1064.25768 
1082.52165 
1098.84984 
1113.53301 
1126.80676 
1138.86395 


1185.65176 
1217.71174 
1241.04485 
1258.78463 
1272.72579 
1283.96982 
1300.98878 
1313.25830 
1322.52269 
1329.76539 


1335.58299 


6.29967 
19.93359 
41.61853 
71.14548 

107.77536 
150.52602 
198.35269 
250.25000 
305.30315 
362.70848 


421.77640 
481.92524 
542.67075 
603.61413 
664.43015 
724.85607 
784.68175 
843.74104 
901.90435 
959.07240 


1015.17085 
1070.14595 
1123.96078 
1176.59214 
1228.02807 
1278.26562 
1327.30917 
1375.16891 
1421.85968 
1467.39994 


62720 10955 
77220 50754 
93880 26333 
12500 74375 
24362 13992 
05527 49635 
32878 63827 
56593 11194 
59738 37326 
35000 00000 


32368 64000 
10754 45816 
65364 43149 
67099 60937 
35099 77590 
11957 76000 
46558 64197 
08705 89636 
32061 92016 
82437 63654 


93247 84000 


63870 95900 
37500 00000 
05861 49151 
99744 15422 
00000 00000 
48184 20411 
94371 30558 
00000 00000 
51688 87156 
25000 00000 


26893 55627 
36128 36364 
27929 25775 
57421 87500 
63289 60000 
43754 31056 
58299 03978 
51653 32896 
93786 73385 
00000 00000 


48575 58614 
76934 57602 
09158 13794 
89338 88081 
33503 89718 
50000 00000 
05359 61070 
25983 20775 
24527 34301 
24218 75000 


ene got em 
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41 
42 
43 
44 
45 
46 
47 
48 
49 
50 


55 
60 
65 
70 
75 
80 
85 
90 
95 
100 


110 
120 
130 
140 
150 
160 
170 
180 
190 
200 


250 
300 
350 
400 
450 
500 


700 
900 
1000 


11 
12 
13 
14 


16 
17 
18 
19 
20 
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Coefficients for Partial Summation of Series—Continued 


A,,(n) 


—1511.81096 
—1555.11614 
— 1597.34045 
— 1638.50997 
— 1678.65152 
—1717.79234 
—1755.95987 
—1793.18151 
— 1829.48448 
— 1864.89568 


— 2029.48358 
—2175.50458 
— 2305.64429 
— 2422.17951 
—2527.02059 
— 2621.76559 
— 2707.75237 
— 2786.10416 
— 2857.76779 
— 2923.54503 


— 3040.07011 
— 3140.04712 
—3226.72629 
— 3302.56973 
—3369.47333 
— 3428.91887 
— 3482.08015 
— 3529.89848 
— 3573.13712 
— 3612.42124 


—3764.93651 
—3869.50853 
— 3945.64787 
—4003.55307 
—4049.06957 
—4085.78689 
—4141.37340 
—4181.45574 
—4211.72539 
—4235.39225 


—4254.40401 


5.17907 
18.43621 
41.05836 
73.11239 

113.91936 
162.42187 
217.42549 
277.74194 
342.26578 


25685 
20836 
08412 
07506 
00000 
61570 
64710 
56936 
33659 
38758 


59023 
12500 
22251 
01945 
71046 
87197 
05561 
00000 
50603 
26601 


23086 
95703 
48747 
16104 
53881 
18030 
64228 
75000 
52882 
64694 


67517 
94386 
09362 
59989 
97094 
38768 
72762 
72446 
76294 
03513 


06764 


08872 
39917 
05027 
36455 
79012 
50000 
46718 
31319 
09298 


15978 
55620 
69187 
36726 
00000 
18721 
78889 
64551 
34024 
40000 


76491 
00000 
61592 
61789 
40000 
87598 
56671 
00000 
22496 
60000 


36282 
12500 
29867 
15303 
60000 
69115 
31147 
00000 
55523 
90000 


07996 
40000 
78951 
80411 
40000 
77630 
97500 
93906 
68478 
60000 


13169 


00610 
69547 
66942 
04849 
34567 
00000 
12227 
75139 
96923 


410.00960 00000 00000 





G0 CO 00 CO GO GO CO 00 CO 00 G0 00 GO 60 CO CO CO CO CO CO GO 00 OO OO 00 00 GO CO CO 00 G0 60 GO 00 00 OO Oo 00 CO 00 00 00 00 60 60 00 60 G0 60 co & 


n 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 


31 
32 
33 
34 
35 
36 
37 
38 
39 
40 


41 
42 
43 
ad 
45 
46 
47 
48 
49 
50 





A, (n) 


480.11474 
551.84849 
624.59437 
697.83950 
771.16137 
844.21501 
916.72132 
988.45671 
1059.24411 
1128.94525 


1197.45405 
1264.69116 
1330.59921 
1395.13900 
1458.28615 
1520.02845 
1580.36354 
1639.29707 
1696.84109 
1753.01280 


1807.83338 
1861.32721 
1913.52101 
1964.44333 
2014.12397 
2062.59359 
2109.88338 
2156.02475 
2201.04914 
2244.98779 


2449.46022 
2631.18696 
2793.38213 
2938.79536 
3069.74844 
3188.19287 
3295.76840 
3393.85582 
3483.62188 
3566.05655 


3712.18554 
3837.65559 
3946.50395 
4041.79455 
4125.89019 
4200.63980 
4267.50953 
4327.67614 
4382.09455 
4431.54741 


74070 


79307 
97129 65006 


18688 37833 
61728 39506 
79968 00000 
65875 63335 
33268 10422 
44642 11340 
79689 38882 
10288 06583 


67567 
21093 
96419 
50920 
60064 
09475 18163 
27443 58135 
20646 08589 
89674 22084 
00000 00000 


51638 
75000 
70076 
20658 
25893 


68884 03656 
00763 42383 
64360 86679 
55667 68517 
33509 45825 
49896 58217 
26923 72969 
56584 36215 
20334 77265 
42272 00000 


07567 22460 
29629 62962 
83334 62128 
02189 56387 
88518 58437 
50000 00000 
32449 16669 
63679 31718 
59619 52758 
77728 00000 


88848 53527 
17695 47326 
00506 44225 
66490 15291 
99969 97532 
27343 75000 
06496 85560 
58449 76206 
34093 57983 
58592 00000 
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Coefficients for Partial Summation of Series—Continued 
m n A,,(n) m n A,,(n) 
8 250 4623.64208 66077 36013 7 55 —1435.20542 13476 13771 
8 300 4755.44125 32133 39918 7 60 —1544.77857 48456 79012 
8 350 4851.44900 77387 00876 7 65 —1642.72273 85062 05653 
8 400 4924.48880 83248 00000 7 70 —1730.64320 00000 00000 
8 450 4981.91638 17457 12963 7 75 —1809.90480 62755 46775 
8 500 5028.25129 74297 95635 7 80 —1881.66054 19464 11133 
8 600 5098.41330 25147 25926 7 85 —1946.88304 53757 79143 
8 700 5149.01707 31826 52370 7 90 —2006.39392 39343 59599 
8 800 5187.23858 37385 75000 7 95 —2060.88921 05422 61924 
8 900 5217.12642 87430 44780 7 100 2110.96075 96896 00000 
8 1000 5241.13788 92728 76238 7 110 —2199.78216 96436 08095 
7 120 —2276.10677 01222-77950 
7 11 — 2.32434 27688 91391 7 130 —2342.36362 83059 88491 
7 12 — 8.82569 23010 97393 7 140 —2400.39978 75000 00000 
7 13 — 20.47422 22035 30323 7 150 —2451.64185 89274 60082 
7 14 — 37.50000 00000 00000 7 160 —2497.20764 05074 59641 
7 15 — 59.64755 88477 36626 7 170 —2537.98449 49786 28544 
7 16 — 86.39309 88311 76757 7 180 —2574.68529 23793 51611 
7 17 — 117.09545 29845 15548 7 190 —2607.88896 80666 75167 
7 18 — 151.09004 10393 62789 7 200 —2638.07036 27369 00000 
7 19 — 187.74223 13549 61766 
7 20 — 226.47432 50000 00000 7 250 —2755.37272 67763 71609 
7 300 —2835.91432 77459 17147 d 

7 21 — 266.77640 60356 65295 7 350 —2894.61302 33856 00000 
7 22 — 308.20785 11549 98331 7 400 —2939.28518 84210 26758 4 
7 23 — 350.39378 76510 47240 7 450 —2974.41821 34306 25530 
7 24 — 393.01911 02832 43312 7 500 —3002.77099 76467 46726 
7 25 — 435.82159 09376 00000 7 600 —3045.71390 96231 66357 
7 26 — 478.58494 40075 21326 7 700 —3076.69358 81384 00000 
7 27 — 521.13230 19392 50302 7 800 —3100.09692 66855 50437 
7 28 — 563.32031 25000 00000 7 900 —3118.39992 63337 66844 
7 29 — 605.03392 90580 70673 
7 30 — 646.18188 75171 46776 7 1000 —3133.10582 56734 31977 
7 31 — 686.69282 35027 79364 | 6 11 .55305 80092 92369 
7 32 — 726.51196 54834 27048 | 6 12 2.18750 00000 00000 
7 33 — 765.59833 42471 88752 6 13 5.21964 71830 56135 
7 34 — 803.92238 08122 51640) 6 14 9.75932 64711 13226 
7 35 — 841.46400 00000 00000 | 6 15 15.76960 00000 00000 
7 36 — 878.21086 35412 96212 6 16 23.12608 33740 23438 
7 37. — 914.15702 35815 42751 6 17 31.66123 64733 16761 
7 38 — 949.30174 41845 75904 | 6 18 41.19341 56378 60082 
7 39 — 983.64852 46062 22385 6 19 51.54427 01561 90718 
7 40 —1017.20428 35937 50000 | 6 20 62.54820 00000 00000 
7 41  —1049.97867 87479 04978 6 21 74.05676 20634 25954 
7 42 —1081.98353 90946 50205 6 22 85.94002 97251 97157 
7 43 —1113.23239 25151 28900 | 6 23 98.08621 69173 04425 
7 44  —1143.74007 26455 01622 6 24 110.40039 06250 00000 
7 45 —1173.52239 23543 72358 | 6 25 122.80277 23776 00000 
7 46 —1202.59587 30014 44941 6 26 135.22692 19685 30348 
7 47 —1230.97752 04371 92854 | 6 27 147.61796 69991 58892 
7 48  —1258.68464 01727 41214 6 28 159.93096 25453 68001 
7 49 —1285.73468 53776 91268 6 29 172.12941 79049 17686 
7 50 —1312.14513 23904 00000 6 30 184.18400 00000 00000 














m 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
6 
§ 


250 
350 
450 
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A,,(n) 
196.07140 79066 45065 


207.77340 
219.27599 
230.5687 1 
241.64404 
252.49691 
263.12425 
273.52468 
283.69815 
293.64575 


303.36947 
312.87201 
322.15667 
331.22719 
340.08767 
348.74247 
357.19616 
365.45341 
373.51901 
381.39779 


418.15736 
450.95050 
480.30905 
506.69735 
530.51306 
552.09389 
571.72594 
589.65167 
606.07704 
621.17759 


647.98371 
671.03743 
691.06399 
708.61598 
724.12092 
737.91423 
750.26251 
761.38011 
771.44130 
780.58909 


816.16379 
840.60878 
858.43361 
872.00420 
882.68005 
891.29754 
904.35278 
913.77345 
920.89153 


53134 
44478 
42727 
64092 
35802 
70025 
08323 
33762 
62500 


39702 
76117 
71429 
78991 
73662 
90551 
07596 
49169 
91578 
54304 


37983 
00000 
71792 
98143 
14784 
39208 
37405 
40740 
65934 
20896 


07590 
43750 
69553 
97483 
71296 
88858 
15381 
93415 
69448 
45944 


99839 
09984 
34170 
72601 
51919 
13687 
57266 
18790 
96405 


91822 
34424 
58785 
34249 
46915 
71157 
13461 
78200 
00000 


45024 
09408 
39398 
97374 
55144 
97216 
15901 
92187 
77238 
00000 


67655 
00000 
66882 
63062 
00000 
98439 
64761 
74074 
56291 
00000 


42449 
00000 
38402 
19153 
00000 
79517 
25483 
63786 
01395 
09000 


54330 
00000 
86880 
23437 
14403 
64007 
00000 
75909 
57838 


926.45915 49364 14814 
930.93306 05482 71897 
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n 

11 
12 
13 
14 
15 
16 
17 
18 
19 
20 


21 
22 
23 
24 
25 
26 
27 
28 
29 
30 


31 
32 
33 
34 
35 
36 
37 
38 
39 
40 


41 
42 
43 
44 
45 
46 
47 
48 
49 
50 


55 
60 
65 
70 
75 
80 
85 
90 
95 
100 


Coefficients for Partial Summation of Series—Continued 





A,,(n) 


.06173 

25117 

61181 
1.16209 
1.90123 
2.81631 
3.88787 
5.09375 
6.41162 
7.82031 


9.30054 
10.83525 
12.40957 
14.01077 
15.62803 
17.25227 
18.87592 
20.49273 
22.09755 
23.68621 


25.25535 
26.80229 
28.32496 
29.82173 
31.29142 
32.73317 
34.14641 
35.53080 
36.88621 
38.21264 


39.51026 
40.77932 
42.02016 
43.23321 
44.41892 
45.57781 
46.71042 
47.81733 
48.89911 
49.95636 


54.89570 
59.31061 
63.26944 
66.83251 
70.05185 
72.97194 
75.63059 


93360 
34825 
72482 
02005 
45679 
94656 
08125 
68798 
66629 
25000 


74154 
34854 
52213 
08212 
20000 
52609 
95975 
43796 
42836 
39917 


29972 
95442 
20804 
63459 
61914 
63133 
52482 
71497 
10502 
64843 


45074 
32831 
75360 
12730 
32671 
48753 
98302 
56978 
67455 
80000 


15874 
92129 
90293 
02635 
60658 
38552 
54008 


99633 
10288 
06838 
11692 
01234 
37208 
08086 
79423 
16569 
00000 


53241 
85399 
16771 
77007 
00000 
43617 
18009 
48361 
86197 
69546 


34498 
15202 
18426 
42708 
67842 
16060 
05055 
26880 
51009 
75000 


77363 
41136 
22183 
48883 
17140 
50937 
11995 
49101 
22113 
00000 


70033 
62962 
69093 
80650 
43621 
85645 
87306 


78.05997 56134 73555 
80.28748 40838 88236 
82.33652 70000 00000 
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2 


110 
120 
130 
140 
150 


170 
180 


200 
250 
350 
400 
450 
500 


700 
800 


3 3 


11 
12 
13 


15 
16 


18 
19 
20 


21 
22 
23 
24 
25 
26 
27 
28 


30 


31 
32 
33 
34 
35 
36 
37 


39 
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Coefficients for Partial Summation of Series—Continued 


A,,(n) 
85.97667 46671 43836 
89.10996 52576 83898 
91.83377 87449 22369 
94.22247 17100 76158 
96.33366 95308 64198 
98.21264 97477 29302 
99.89544 02988 96711 
101.41104 24975 44413 
102.78304 95989 64891 
104.03084 19062 50000 


108.88638 64197 12000 
112.22555 24084 36215 
114.66175 30495 62682 
116.51724 88093 87207 
117.97739 09714 30507 
119.15629 02246 08000 
120.94275 99181 13426 
122.23222 49461 38513 
123.20671 31642 56095 
123.96905 07440 93888 


124.58170 56530 97000 


.00231 20852 17500 
00960 21947 87380 
.02376 06252 90953 
04569 52460 28441 
07551 47325 10288 
.11279 29687 50000 
-15679 72317 34479 
-20665 32314 97005 
.26145 30270 99227 
.32032 00000 00000 


38244 43453 61032 
44709 94789 34115 
51364 66873 92137 
.58153 29218 10700 
-65028 48921 60000 
-71950 14345 91672 
-78884 53416 30860 
85803 53424 16850 
.92683 86032 22821 
-99506 39231 82442 


1.06255 56830 78909 
1.12918 85375 97656 
1.19486 28049 94586 
1.25950 04907 08074 
1.32304 18754 09055 
1.38544 25985 19872 
1.44667 11723 35430 
1.50670 68677 06442 
1.56553 79186 90180 
1.62316 00000 00000 





-_ PEHLALALALAHLE PEELE EE EE PEP EEE EE ES -p-_DLEE Ee eS S 


n 

41 
42 
43 
44 
45 
46 
47 
48 
49 
50 


55 
60 
65 
70 
75 
80 
85 
90 
95 
100 


110 
120 
130 
140 
150 
160 
170 
180 
190 
200 


250 
300 
350 
400 
450 
500 
600 
700 
800 
900 


1000 


A,,(n) 
1.67957 49371 42942 
1.73478 96146 54486 
1.78881 50527 58094 
1.84166 56270 93845 
1.89335 84099 08155 
1.94391 26143 25570 
1.99334 91260 90177 
2.04169 01095 25035 
2.08895 86764 66958 
2.13517 86086 40000 


2.35140 17642 06820 
2.54504 69135 80248 
2.71896 82698 94251 
2.87571 51567 79913 
3.01750 50773 63183 
3.14624 29687 50000 
3.26355 66750 23902 
3.37083 54864 60396 
3.46926 63012 67056 
3.55986 54873 60000 


3.72094 07034 81281 
3.85971 10836 76269 
3.98043 55316 31768 
4.08637 31730 82645 
4.18005 42836 72757 
4.26346 98852 53906 
4.33820 62246 81532 
4.40554 15310 44839 
4.46651 66258 98876 
4.52198 71590 40000 


4.73798 01192 37264 
4.88664 43058 07626 
4.99516 96792 10730 
5.07786 09345 10000 
5.14295 36638 91081 
5.19552 17538 96715 
5.27520 39203 87161 
5.33273 43969 14505 
5.37622 08939 50156 
5.41024 55059 97484 


5.43759 28570 10013 
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model 402 accounting machine, and the sum checks 


(3) XAn(m) = 1 


were totalled and listed concurrently. Inspecting the sum checks and spot-check- 
ing by hand-computation suggest that the results A,,(m) are accurate to 14 
decimals. 

The A,,(m) were rounded on cards to 15 decimals, using an International 
Business Machines Corporation model 075 sorter for machine inspection of the 
16th decimal, and an International Business Machines Corporation model 513 
reproducer for subsequent insertion of the corrected 15th decimal. It is believed 
that the error in A,,(m) is < 1 X 107. 


R. B. HorGAN 
University of California 
Los Angeles, California 


1. HerBert E. Sauzer, “Formulas for the partial summation of series’”’ [see p. 149 this issue ]. 


Polynomial Approximations to Some Modified Bessel Functions 


The following approximations were determined by equating tabular values of 
the functions to polynomials of degree m at nm + 1 points and solving the equations 
for the m + 1 coefficients. In the region from 0 to x,, the fit points were spaced 
at equal intervals in x*; in the region from x, to ©, they were spaced at equal 
intervals in 1/x. Both m and x, were determined by trial and error with the 
following requirements. 


(1) To yield 7 or 8 accurate significant figures over the full range of positive x. 

(2) To make the polynomial degrees approximately the same in the two 
regions 0 to x, and x, to ~. 

(3) To keep the magnitudes of the polynomial coefficients, the definitions 
of the polynomials, and the regions of validity balanced to the extent that a given 
number of significant figures in the coefficients would give the same number of 
significant figures in the polynomials over the entire region of validity. 


In order to determine the maximum errors, the polynomials were evaluated 
and compared (directly or after appropriate multiplications and additions) with 
values given in [1]. Differences between calculated and tabular values were 
plotted for 

x = .O[.1]5[.2]7.6[.4 ]10[1 ]20. 


To determine maximum error estimates in formulas (3) and (4), it was 
necessary to compute J» and J; at higher values of x from asymptotic expansions. 
Local maxima in these two formulas occur approximately at midpoints between 
fit points, decreasing gradually as x is decreased. The lowest is about .1 X 10-7 
at x = 4.1 and the highest, about 2 X 10-7 at x = 60. 

These approximations may be considered a sequel to a similar set for Jo, Yo, 
J, and Y;, which was published earlier [2]. 
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(1) — 3.75 <x < 3.75 || max ~ 3.0 X 10-* 


Io(x) ~ 1 + 3.5156 229 (x/3.75)? + 3.0899 424 (x/3.75)* 
+ 1.2067 492 (x/3.75)* + .2659 732 (x/3.75)*® 
+ .0360 768 (x/3.75)" + .0045 813 (x/3.75)" 


(2) — 3.75 <x < 3.75 | max ~ 1.0 X 10-* 


Ii (x)/x = 5 + .8789 0594 (x/3.75)? + .5149 8869 (x/3.75)* 
+ .1508 4934 (x/3.75)* + .0265 8733 (x/3.75)*® 
+ .0030 1532 (x/3.75)" + .0003 2411 (x/3.75)" 


(3) 3.75<x< @ |e] max ~ 2 X 10-7 


To(x)xte-= ~ 3989 42280 + .0132 85917 (3.75/x) 
0022 53187 (3.75/x)? — .0015 75649 (3.75/x)* 
0091 62808 (3.75/x)* — .0205 77063 (3.75/x)® 
0263 55372 (3.75/x)* — .0164 76329 (3.75/x)? 
.0039 23767 (3.75/x)* 


++++ 


(4) 3.75<x*%< @ le| max ~ 2 X 10-7 


Ii(x)x¥e-= ~ 3989 42280 — .0398 80242 (3.75/x) 
— .0036 20183 (3.75/x)? + .0016 38014 (3.75/x)* 
— .0103 15550 (3.75/x)* + .0228 29673 (3.75/x)* 
— .0289 53121 (3.75/x)* + .0178 76535 (3.75/x)? 
— .0042 00587 (3.75/x)* 


(5)0<x<2 \e| max ~ 7.0 X 10-8 


Ko(x) + In (.5x) Io(x) ~ — .5772 1566 + 4227 8420 (x/2)? 
+ .2306 9756 (x/2)* + .0348 8590 (x/2)* 
+ .0026 2698 (x/2)* + .0001 0750 (x/2)" 
+ .0000 0740 (x/2)" 


(6) 0<zx<2 le| max ~ 6.0 X 10-8 
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{K.i(x) — In (.5x) Ii(x)}x~ 1 + .1544 3144 (x/2)? 
— .6727 8579 (x/2)* — .1815 6897 (x/2)* 
— .0191 9402 (x/2)* — .0011 0404 (x/2)" 


— .0000 4686 (x/2)" 


(7)2<x<o@ le| max ~ 1.5 X 10-7 


Ko(x)xtes~ 1.2533 1414 — .0783 2358 (2/x) 
+ .0218 9568 (2/x)* — .0106 2446 (2/x)* 
+ .0058 7872 (2/x)* — .0025 1540 (2/x)* 
+ 0005 3208 (2/x)* 


(83) 2<x<o@ \e| max ~ 1.5 X 10-7 


Ki(x)xte?~ 1.2533 1414 + .2349 8619 (2/x) 
— .0365 5620 (2/x)? + .0150 4268 (2/x)* 
— .0078 0353 (2/x)* + .0032 5614 (2/x)* 
— .0006 8245 (2/x)* 
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Iculated F — tabular F, 
In these formulas, « = oe — where F = Io, J1, Ko, or 
tabular F 
K,. [e given in the approximations to Jo, Yo, Ji, and V1 previously mentioned was 


the absolute error in the functions represented by polynomials. ] 


E. E. ALLEN 





Shell Development Company 
Exploration and Production Research Division 
Houston, Texas 


1. BritisH ASSOCIATION FOR THE ADVANCEMENT OF SCIENCE, Mathematical Tables, v. X, 
Bessel Functions, Cambridge Univ. Press, 1952. 
2. E. E. ALLEN, “Analytical approximations,” MTAC, v. 8, 1954, p. 240-241. 


Radix Tables for sin x and cos x, x = a-10* degrees, 
a= 1(1)9,k = —3(1)1 


These tables permit calculation of sin x or cos x for any argument x expressed 
to 3D by means of summation formulas listed in Table I and requiring at most 
64 multiplications and 15 additions or subtractions. 


Maximum number of 


No. of decimals Maximum number of additions or 
in x multiplications subtractions 

3 64 15 

2 24 7 

1 8 i 3 

0 2 1 


The table was computed on the SWAC machine using double precision rou- 
tines. (Each SWAC number has a sign plus 36 binary digits.) Values of 1/n! for 
m = 2(1)20 were taken from Peters and Stein [1 ]. 

For x = .001(.001).01(.01).1(.1)1(1)10(10)50°, expressed in radians, double 
precision, the routine computed sin x and cos x by 


. Ce. #. - Od 
sins = %— 3 a +++ = oe 
and 
Rr Sgioigt x 
RES = EA TIRT Het F a7 


For each x, as a term — was computed, it was accumulated into the sin x or the 
m! 


cos x cell (alternating) with the proper sign attached. The limit of 20 was chosen 
for m in the formulas above, since there would be no contribution from succeeding 








st 


or 


le 
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terms in the power series for the functions. The values for sin x and cos x where 
x = 60(10)80° were lifted from the cos x and sin x for x = 10(10)30°. 

The values were converted by SWAC to decimal form, tabulated on an 
International Business Machines Corporation model 402 accounting machine, 
and rounded to 20D by hand. Accuracy to 0.6 X 10-* is believed assured. 


The tables were computed in connection with the review of AMS 40 in this 
issue, p. 169. 
R. B. HorGAN 
University of California 
Los Angeles, California 


1. J. Peters & J. Sremn, Appendix to Zehnstellige Logarithmentafel. Band I. Zehnstellige 
Logarithmen der Zahlen von 1 bis 100,000 nebst einem Anhang mathematischer Tafeln, Berlin, 
Reichsamt fiir Landesaufnahme, 1922. A reprint is being issued by F. Ungar, New York. 


TABLE I. Summation Formulas for Sine and Cosine 


sin (a + 8) = sin acos 8 + cosasin B 
cos (a + 8) = cosacos 8 — sina sin B 


sin (a + 8 + y) = sinacos #8 cos y — sinasin 8 sin y + cos acos 8 sin y 
+ cos a sin 6 cos y 


cos (a + 8 + 7) = cosacos # cos y — cosasin 8 sin y — sin asin 6 cos y 
— sinacos 8 sin y 


sin (a + 6 + y + 8) = sina cos 6 cos 7 cos 6 — sin asin @ sin y cos 6 
— sin asin 8 cos y sin 6 — sin acos 8 sin y sin 6 
+ cos a sin 8 cos y cos 6 + cos a cos B sin y cos 6 
+ cos a cos 6 cos y sin 6 — cosa sin # sin y sin 6 


cos (a + 8 + y + 8) = cosacos# cos y cos § — cos asin B sin y cos 6 
— cos asin 8 cos y sin é — cosacos@ sin y sin é 
+ sin a sin 8 sin y sin 6 — sin a sin 6 cos y cos 6 
— sin acos 8 sin y cos 6 — sin acos § cos y sin 6 


sin (a+68++7+6+6) =sinasin£ sin y sin é sin e — sin a cos 8 cos y sin 4 sine 
— sin acos 8 sin y cos é sin e — sin asin Bcos7 cosésine 
— sin acos 8 sin y sin 6 cos e — sin asin8cos7ysinécose 
— sinasin # sin y cos dcose + sinacos# cosy cosécose 
— cosacos sin y sin é sin e — cosasinScosysinésine 
— cos asin # sin y cos sine + cosacosfcos7cosésine 
— cos asin # sin 7 sin cose + cosacos#cos 7 sin 6 cose 
+ cos a cos Ssin y cosécose + cosasin 8 cos 7 cos 6 cos € 


cos (a + 8 + 7 +6 +) = cosacos # cos 7 cos 6 cos e — cos asin # sin y cos 6 cos € 
— cos asin B cos y sin 6cose — cosacos§ sin 7 sind cose 
— cos asin 6 cos y cos dsine — cosacos#sin 7 cosé sine 
— cos a cos 6 cos y sin é sin e + cosasin# sin ysinésine 
— sin asin B cos y cos 6cose — sinacos#sin y cosé cose 
— sin acos cos y sin 6 cose + sinasin#sinysinédcose 
— sin a cos 6 cos y cos 6 sin e + sinasin#sin y cosésine 
+ sin asin 6 cos y sin 6 sin e + sin acos@ sinysinésine 
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sin x 


00001 74532 
-00003 49065 
00005 23598 
-00006 98131 
.00008 72664 
.00010 47197 
.00012 21730 
.00013 96263 
-00015 70796 


.00017 45329 
-00034 90658 
00052 35987 
-00069 81316 
.00087 26645 
.00104 71973 
00122 17301 
.00139 62629 
00157 07956 


.00174 53283 
.00349 06514 
00523 59638 
.00698 12602 
.00872 65354 
.01047 17841 
.01221 70008 
.01396 21803 
.01570 73173 


.01745 24064 
03489 94967 
.05233 59562 
.06975 64737 
08715 57427 
-10452 84632 
.12186 93434 
13917 31009 
-15643 44650 


.17364 81776 
.34202 01433 
-50000 00000 
.64278 76096 
-76604 44431 
-86602 54037 
-93969 26207 
-98480 77530 
1.00000 00000 00000 00000 


92519 05720 
85032 79782 
77535 90529 
70023 06303 
62488 95446 
54928 26301 
47335 67209 
39705 86514 
32033 52557 


24313 33680 
43310 09671 
51673 70300 
44087 57925 
15235 14954 
59799 83861 
72465 07198 
47914 27617 
80830 87881 


65898 30884 
15223 73227 
31419 58009 
97961 55247 
98373 93496 
16245 79346 
35247 16888 
39145 27162 
11820 67575 


37283 51282 
02500 97165 
42943 83272 
44125 30078 
47658 17356 
67653 47140 
05147 48111 
60065 44411 
40230 86901 


66930 34885 
25668 73304 
00000 00000 
86539 32632 
18978 03520 
84438 64676 
85908 38405 
12208 05937 





TABLE II. Radix Tables for Sine x and Cosine x in Degrees 


cos xX 


-99999 99998 
-99999 99993 
-99999 99986 
.99999 99975 
-99999 99961 
-99999 99945 
.99999 99925 
-99999 99902 
.99999 99876 


99999 99847 
-99999 99390 
-99999 98629 
99999 97563 
-99999 96192 
-99999 94516 
-99999 92536 
-99999 90252 
.99999 87662 


-99999 84769 
.99999 39076 
-99998 62922 
-99997 56307 
-99996 19230 
99994 51693 
-99992 53696 
-99990 25240 
-99987 66324 


-99984 76951 
-99939 08270 
-99862 95347 
.99756 40502 
.99619 46980 
-99452 18953 
-99254 61516 
-99026 80687 
.98768 83405 


-98480 77530 
-93969 26207 
-86602 54037 
-76604 44431 
-64278 76096 
.50000 00000 
.34202 01433 
17364 81776 
-00000 00000 


47691 29011 
90765 16049 
29221 61127 
63060 64270 
92282 25508 
16886 44885 
36873 22451 
52242 58266 
62994 52401 


69129 04933 
76516 66127 
22164 22770 
06074 06842 
28249 43114 
88694 49148 
87414 35299 
24415 04715 
99703 53332 


13287 69880 
57790 38148 
47426 79269 
05394 79311 
64171 28874 
€5512 13198 
60451 99446 
09304 21155 
81660 59864 


56391 23916 
19095 73001 
54573 87378 
59824 24761 
91745 53229 
68273 33692 
41322 03498 
41570 31508 
95137 72619 


12208 05937 
85908 38405 
84438 64676 
18978 03520 
86539 32632 
00000 00000 
25668 73304 
66930 34885 
00000 00000 
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Note on Predictor-Corrector Formulas 


The various predictor-corrector formulas for solving differential equations 
have considerable appeal because of their local control of truncation error and 
machine error. Of course, there is an expense for this control in terms of additional 
computing, so that a different method might be chosen if the accuracy does not 
need to be examined at each step. However, predictor-corrector formulas are 
quite generally used [1], [2]. 

There seems to be a natural tendency to repeat the “‘correction’’ formula, 
just like you take two aspirin instead of one. In fact the convergence of the 
iterative procedure defined by a “‘correction”’ formula has been studied, but evi- 
dently without inquiring whether the successive steps in the iteration actually 
afford an improvement. The purpose of this note is to point out that repeating a 
“correction” formula may worsen the result rather than improve it. : 

Consider y’ = f(x, y), initial value (xo, yo), predictor y?n41 = Ya—-1 + 2hf (Xa, Yn), 
and corrector y°ng1 = Yn + (h/2)[f (Xm, Yn) + f(Xn+1, Y?n41) ]- 

This is a commonly used method with truncation errors E? = (h*/3)y’”, 
E* = — (h®/12)y’”; h is the increment in x, and y’” denotes an appropriate mean 
value of y’”’ (x). If y’”’ (x) > 0 near the point in question and y,,; is the theoretical 
value we are seeking, then the signs of E? and E* show that y?asi — Yasi < 0 
while ati — Yori > 0; also Wnt > Vrn+i- 

Now suppose 0f/dy > 0 near the point in question: then applying the “‘cor- 
rector’”’ to the “corrected’’ point gives 


Yngt = Yn + (h/2) Cf (Kn, Yn) + f(Xnt1, Yuet) ] > nga > Yar» 


so that y,1; is farther than y*,4: from yn4:. It is also easy to give a geometric 
illustration for this result. 

Criteria could be given indicating when repetition of the ‘‘corrector’”’ formula 
will worsen and when it will improve y*,;:. The situation is similar with other 
predictor-corrector formulas, and improvement is only a 50-50 proposition in 
that it hinges on the sign of a derivative. However, it seems much more practical 
to never repeat a “corrector,” and if |-y°n41 — y?n41| is too large to accept, then 
a smaller value of h should be used. 


D. D. WALL 
IBM Corporation 
Los Angeles, California 
1. F. B. HiLpEBRAND, Introduction to Numerical Analysis, McGraw-Hill, New York, 1956. 


2. W. E. Mine, Numerical Solution of Differential Equations, John Wiley and Sons, New 
York, 1953. 


The Order of an Iteration Formula 


The order of an iteration formula indicates the rate of convergence of the 
iteration [1]. For example, Newton’s formula x4: = x, — f(x.)/f’(x,), for 
solving f(x) = 0, is of order 2 since the number of correct decimal places approxi- 
Xnz1 — a = (xq, — a)*f’’(a)/2f’ (a) + 0, where @ consists of terms of degree three 
and higher in (x, — a) and when f(a) = 0. More generally, for the iteration 
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formula xn: = $(Xo, ---, Xn) which converges to a, let p, = — log|x, —- a| so 
that p, is the number of correct decimal places in x,; then we define the order p 
of the iteration formula to be p = lim Pn41/p, asn —> ©. 

If ¢ depends on x, only, then the order is an integer. In fact the order in this 
case is the index p of the smallest derivative of @ such that ¢ (a) ¥ 0, as is seen 
by the power series expansion of ¢(x,) in powers of (x, — a). 

We now show that a particular iteration formula has order (1 + V5)/ 2 = 1.6+, 
which is interesting because this number is not an integer and also because it is 
large enough to indicate that the formula is quite useful. The formula in question 
is a common one for solving y = f(x) = 0 by repeated linear interpolation: x41 
is the point where the straight line through (x,_1, yn_1) and (xn, yn) crosses the x 
axis. Algebraically this gives %n41 = (%n—-1¥n — XnYn—1)/ (Yn — Yn—1). The fact that 
the formula does not involve f’ (x) will often be a considerable advantage [2]. 

To investigate the order of the linear interpolation formula we expand 
y = f(x) in powers of (x — a) and with a little algebraic simplification obtain 
Kapa — @ = (Xa-1 — a) (Xn — a) f’’(a)/2f’ (a) + 6, where @ consists of terms of 
degree three and higher. From this equation it follows that Pay: = pn + Pas 
+ qn, where g, is bounded, so on dividing by p, and letting n> © we get 
p = 1 + 1/p, whose only positive root is p = (1 + V5)/2. This is the desired 
result. 

It may also be of interest to note that the interpolation formula can be 
written as Xn41 = (Xn — Yn) + (yn — Yn—1)/(%n — Xn—1) ], which shows a close 
relation to Newton’s formula. Although (yn — yn—1)/(*n — Xn—1) approaches f’ (a), 
the ‘‘one sidedness”’ of this difference quotient keeps the order of the interpolation 
formula less than the order of Newton’s formula. 


- D. D. WALL 
IBM Corporation 


Los Angeles, California 


1. D. R. Hartree, “Notes on iterative processes,” Cambridge Phil. Soc., Proc., v. 45, 1949, 
p. 230-236; see also E. SCHRODER, “Ober unendlich viele Algorithmen ziir Auflésung der Glei- 
chungen,” Math. Ann., v. 2, 1870, p. 317-363. 

2. A. S. HouSEHOLDER, Principles of Numerical Analysis, McGraw-Hill Book Co., Inc., New 
York, 1953, Sec. 3.3, p. 118-132. 
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48[C, F].—F. J. Duarte, ‘Tablas Logaritmicas de Factoriales Primarias, desde 
2 hasta 1007, con 33 Decimales,’’ Acad. Ci. Fis., Mat. y Nat., Bol., Caracas, 
v. 15, 1953. 


This booklet gives, for each prime p = 2, 3, ---, 10007, the value to 33D 
of logio (p) where (p) = 2-3-5---p; the last digit is asterisked when it has 
been rounded up. There is a short table of log. (p)/y, for y = 50, 100(100) 
1000(1000)10000, where # is the greatest prime less than y, exhibiting the ap- 
proach of this ratio to unity, a result which is equivalent to the prime number 
theorem. 

This table is, presumably, based on the 36D table [1] by the same author 
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which gives logio m for m = 1(1)1000(primes) 10007. It is apparently reproduced 
photographically from typescript. 
5p T: 


1. F. J. Duarte, Nouvelles Tables Logarithmiques a 36 décimales, Gauthier-Villars, Paris, 1933. 
[MTAC, v. 1, 1943-45, p. 3.] 


49[D ].—NBS Applied Mathematics Series, No. 40, Table of Secants and Cosecants 
to Nine Significant Figures at Hundredths of a Degree, U.S. Govt. Printing 
Office, Washington, D. C., 1954, vi + 46 p., 26 cm. Price $0.35. 

This table lists sec x and csc x = 0(.1)90°, 9S. 

The table was computed by inverting values in [1]. Differencing and other 
tests were made, and the introduction states a belief that entries are correct within 
six-tenths of a unit in the last place given. 

The introduction provides information concerning interpolation and the use of 
a truncated power series in the region where interpolation is not feasible. 

The table is printed photographically, presumably from masters prepared 
from a card-controlled typewriter (although there is no announcement in the 
preface or in the introduction as to how the printing was done). 

A spot check of 200 values for each of the functions sec x and csc x was under- 
taken on the SWAC computer. 200 random arguments were chosen from [2] for 
each function. The values of sin x and cos x were computed at these 200 values 
from radix tables prepared especially for the purpose—that is, sin x and cos x 
were computed for x = .01(.01).1(.1)1(1)10(10)90°, 21D [3]. 

From these tables, which were stored in the SWAC, the values of sin x and 
cos x for any argument in the range is computable as the function of the sum of 
four arguments in the table. For the arguments chosen either sin x or cos x (as 
appropriate to the function being tested) was calculated and the product of this 
value with the value listed in the table was compared with 1. No divergence was 
found which corresponds to an error as great as six-tenths in the last digit listed. 


RutH B. HorGan 


om s. 

University of California 
Los Angeles, California 

1. NBS Applied Mathematics Series 5, Tables of Sines and Cosines to Fifteen Decimal Places 
at Hundredths of a Degree, U. S. Govt. Printing Office, Washington, D. C., 1954. 

2. The RAND Corporation, A Million Random Digits with 100,000 Normal Deviates, Free 
Press, Glencoe, Illinois, 1955. 

3. R. B. HorGan, “Radix tables for sin x and cos x, x = a-10* degrees, a = 1(1)9,k = —3(1)1” 
[see p. 164 this issue ]. 


50(H, Z ].—Proceedings of the Second Symposium in Linear Programming, held in 
Washington, D. C., January 27-29, 1955, edited by H. A. Antosiewicz. 
National Bureau of Standards, and Directorate of Management Analysis 
DCS/Comptroller, U.S.A.F. vi + 685 p. 

These 2 volumes contain 27 papers in full, and the abstracts of 6 others, pre- 
sented at the linear programming symposium sponsored by the Office of Scientific 
Research of the Air Research and Development Command. The papers are 
grouped into five parts: Applications, Economic Theory, Computation, Theory 
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of Linear Inequalities, and Developments in Linear Programming. Those papers 
most directly concerned with computational methods are described briefly here. 

L. Gainen’s paper, ‘Linear programming in bid evaluations,”’ outlines the 
method of computing cost contract awards on S.E.A.C. The problem is essentially 
the Hitchcock distribution problem, but with extra constraints such as ‘Either 
> Xy = Oor Y Xi > Aj.” This type of constraint involves solving the problem 
j r) 


separately under each condition, and selecting the cheaper of the two solutions. 
Many similar complications are discussed. 

S. Johnson and G. B. Dantzig’s paper, ‘‘A production smoothing problem,”’ 
presents a graphical method of minimizing the cost of production of a single item 
over a given number of time periods to satisfy known future requirements, when 
the cost per unit for production, storage, and change in production rates are 
known functions of time. A good first approximation to the solution is obtained 
by considering the original problem and its dual simultaneously. Often this is 
clearly the required solution, if not, the latter is obtained with a few iterations of 
the simplex method which can easily be performed graphically. 

R. Bellman’s paper, “Dynamic programming and multi-stage decision proc- 
esses of stochastic type,”’ consists of a brief introduction to the author’s theory of 
dynamic programming, in which sequential decision problems are formulated as 
functional equations. 

A. J. Hoffman’s paper, ‘How to solve a linear programming problem,” consists 
of a general introduction to the computational aspects of linear programming. 
Machine methods are discussed in the light of experience on S.E.A.C., the IBM 
701 at RAND, and other machines. The author suggests that some version of 
the simplex method should always be used, but that special devices should be 
incorporated to take advantage of special features. Thus a special routine should 
be used for the distribution problem. 

C. Tompkins’s paper, “Projection methods in calculation,”’ discusses a general 
relaxation method of solving simultaneous polynomial equations or inequalities. 
Methods of overrelaxing and acceleration are indicated. It is claimed that with 
enough acceleration one can solve an extensive system in an acceptable time. 

A. Vazsonyi’s paper, ‘Optimizing a function of additively separated variables 
subject to a simple restriction,’ discusses a simple technique for maximizing 
>X fi(x:), where the f; are general functions and the x; are nonnegative variables 
subject to a single linear inequality. 

R. M. Thrall’s paper, ‘Some results in non-linear programming,”’ is chiefly 
concerned with the special case of the problem considered in the previous paper 
when the f;(x,) are all convex functions. 

S. I. Gass’s paper, “‘A first feasible solution to the linear programming prob- 
lem,”’ shows how such a solution can be found in terms of a single artificial variable 
if the constraints have been solved for m of the variables—preferably those which 
are expected to have positive values in the final solution. The method is in effect 
the dual to that considered in [1]. 

H. Markowitz’s paper, ‘Concepts and computing procedures for certain Xj; 
programming problems,” is primarily concerned with problems of the follow- 
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ing form— 


Xu ta = Mi, 
i 


LD aX i 


v9. 


Minimize >> CXi;, where ai, Ci; > 0, the variables X;; and ¢; are restricted 
to nonnegative values, and @ is a parameter. Properties of the basic optimal 
feasible solutions of this type of problem are established, and an algorithm, based 
on the simplex method, for finding optimal solutions as @ varies is described. 
Extensions to more general problems are discussed. 

E. D. Schell’s paper, “Distribution of a product by several properties,”’ 
discusses the distribution problem in 3 dimensions, i.e., the minimization of 
X Pints for nonnegative x; subject to constraints of the form >> xi = uy or 

i 


xX >} xis = a;. Methods of finding feasible solutions, and iterative methods, 
3 tf 

based on the simplex method, of proceeding to the optimal solution are discussed. 
It is shown that the optimal solution may involve fractional values of the variables 
even when the coefficients in the constraints are all integers. 

I. Heller’s paper, “On the traveling salesman’s problem,” provides theoretical 
background to some semi-empirical methods of solving this problem used suc- 
cessfully at the RAND Corporation. 

In G. B. Dantzig’s concluding paper, ‘‘Developments in linear programming,” 
a method of solving the distribution problem when the variables have fixed upper 
bounds is indicated. It is shown that this method can be used to solve the problem 
when the requirements at each destination are unknown but have a known 
probability distribution. 

Possible developments in methods and fields of application are reviewed in 
more general terms. 


E. M. L. BEALE 
Admiralty Research Lab. 
Teddington, Middlesex 
England 


1. E. M. L. Beate, “An alternative method for linear programming,’’ Camb. Phil. Soc., Proc., 
v. 50, 1954, p. 513-523. 


51[1].—J. Kuntzmann, ‘“‘Formules de Dérivation Approchée au Moyen de Points 
Equidistants,”” Report No. 1.373/1, Société d’Electronique et d’ Automatisme, 
138 Boulevard de Verdun, Courbevoie (Seine), France, 1954, 30.5 cm. 


These tables give the coefficients A; »., in the finite difference approximation 
i=n 
of the derivatives of f(x); f(x.) = & Ax »:ef (xi) for p = 1(1)m, g = O(1)m; 


n = 2(1)10. The tables are extensions of the earlier work [1] of Bickley. The 
coefficients are given as rational fractions with a common denominator. The 
coefficients for the first derivative were computed directly and the balance were 
obtained recursively from the relation M,-M, = M,,, where M, is the matrix 
whose elements are A; »;,- Checks were made by applying the results to the func- 
tion x". 
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Values at intermediate points, for the first derivative [2], in the case of 
n = 4(1)7, have been given by H. E. Salzer. 


MILTON ABRAMOWITZ 
National Bureau of Standards 


“re - ¢ 


. W. G. Bicxiey, “Formulae for numerical differentiation,” Mathematical Gazette, v. 25, 
14 a 19-27. 


E. Satzer, Table of Coefficients for obtaining the First Derivative without Differences, 
National Bureau of Standards, Applied Mathematics Series 2, 1948. 


52[K ].—W. ALLEN Watuis & Harry V. Roserts, Statistics: A New Approach, 
The Free Press, Glencoe, Illinois, 1956, xxxviii + 646 p., 23.4 cm. Price $6.00. 
Includes the first 10,000 digits from the “RAND Random Digits’’ [1]. 

te Sen, Sen 


1. The RAND Corporation, One Million Random Digits and 100,000 Normal Deviates, The 
Free Press, 1955, Review 11, MTAC, v. 10, 1956, p. 39-43. 


53[L].—J. A. Stratton, P. M. Morsg, L. J. Cuu, J. D. C. Lirrie, & F. J. 
CorBato, Spherical Wave Functions, Including Tables of Separation Constants 
and Coefficients, The Technology Press, Mass. Inst. of Tech., and John Wiley 
& Sons, Inc., New York, 1956, xiii + 613 p. Price $12.50. 

Angular and radial spheroidal wave functions are solutions of second order 
ordinary differential equations which arise when the wave equation is solved by 
using the separation of variables in spheroidal coordinates. The set of tables 
contained in this book gives the coefficients of the expansions of the radial func- 
tion of the first kind in spherical Bessel functions and the angular function of the 
first kind in associated Legendre functions. Tables for both the prolate and oblate 
cases are given. In these expansions two integers, m and /, occur with / > m and 
a parameter h or g. The tables cover the ranges 


m=0(1)8 1 = m(1)8 
g, h = 0(.1)1(.2)8. 
The table entry 


+0.6278567 | — 02 means +0.6278567 x 10~. 


The authors state that: 

“The accuracy of the tables has been checked by the hand calculation of a 
number of cases. For the most part, all the significant figures published are good 
except for an occasional error of one in the last place. When the coefficient itself 
gets large (~10*), the last place may be weak. At no time is it thought that the 
numbers are off by more than 5 in the last place.” 

In addition to the actual tables, this book contains a foreword by P. M. Morse 
and a reprinted version of the paper by L. J. Chu and J. A. Stratton [1]. In this 
paper, a variety of properties of the spheroidal wave functions are derived and 
discussed. 

The book also contains a section entitled ‘“‘Introduction to the Tables’ by 
John D. C. Little and Fernando J. Corbato. In this section the number notation 
of the tables is explained, the accuracy of the numbers is discussed, the sources 
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of auxiliary tables are listed, and a brief description is given of the programming 
of the table calculations for the Whirlwind I computer. 

Morse states in his introduction that “. . . the present tables were computed, 
tabulated and printed automatically ; no human intervention entered between the 
introduction of the program and the typed page, ready for offset photography.” 

One can thus expect a considerable reduction in errors, in part because of the 
inherently greater accuracy of the machine over human operations with desk 
calculators and in part because of the fewer transcriptions between calculated 
results and final printed table. 

A. H. T. 


1. L. J. Cau & J. A. Srratton, ae and spheroidal wave functions,” J. Math. Phys., 
Mass. Inst. Tech., v. 20, 1941, p. 259- 


54(L ].—Cecix Hastincs, Jr., “Approximations for calculating Fresnel Integrals,” 
Cecil Hastings’ “Approximation Newsletter,” April 12, 1956, Note-10. 
The functions (7°) and A(7*) are approximated to about 3D, 0 < T < o, 
where 
T iA iz mi edt 
(x) + tA) = - Vost 


Multiplication of these functions by sines and cosines yield the Fresnel In- 
tegrals in accordance with classical formulas [1]. Note that modern notation 
for C(x) and S(x) as used in [2] and [3] differs from that used here (see MTAC, 
v. 9, 1955, p. 124). 

The approximations are 


I(T?) ~ (2 + 3.305T + 2.2237? + 3.3887), 
A(T?) ~ (1 + 0.7397) (2 + 1.430T + 1.97672). 


For tables of the Fresnel Integrals see [2] and [3]. 
ca 2 
EvuGENE JAHNKE & FRITz are Tables of Functions with Formulae and Curves, Dover 
Publications New York, 1945, 
2. A. VAN WIJNGAARDEN ew W. L. ScHEEN, “Table of Fresnel Integrals," Akademie van 


Wetenschappen, Amsterdam, Afd. Natuurkunde, Verhandelingen, eerste sectie, v. 19, no. 4, 1949 
[MTAC, v. 4, 1950, p. 155-156]. 


3, AKADEMIIA NAUK SSSR. Tabliisy Integralov Frenelya, Moscow, 1953 [MTAC, v. 9, 
1955, p. 75-76]. 


55[L].—E. A. Taytor & J. M. KenneEpy, Tables of Complete Elliptic Integrals, 
Atomic Energy of Canada Limited, Chalk River, Ontario, A.E.C.L. Report 
No. 296, 1956, iii + 10 p., 27.5 cm. Price 50¢. 

These tables give values of the well known complete elliptic integrals K and 

E to 6D for k* = 0(.001)1. There are no differences. The bulk of the values were 

obtained by interpolation in tables to 7D by Milne-Thomson, but some values 

at each end of the range were computed independently. The authors state that 
the errors should never exceed a unit in the last place. 

Values to at least as many decimals for the same arguments have been pub- 
lished by several other authors, though in works which may possibly rank as not 
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highly accessible. In fact, values to 10 or more decimals were printed in two tables 
published by Hayashi in 1930 and 1933, and the few errors which exceed a unit 
or two in the last place are all known, so that a table to 6D and entirely free from 
rounding errors could be produced by anyone with access to the Hayashi tables 
and the corrections, without computing any values. For details, see MTAC, v. 3, 
1948, p. 232, 263. 

Sampling by comparison of one-fifth of the values with those of Hayashi sug- 
gests that the number of slight rounding errors averages approximately two per 
column of 50 values. Evidently these will not affect the utility of the table. 

A. F. 


56[L ].—Lucy J. Starter, “A short table of the Laguerre polynomials,” Inst. 
Elec. Engineers, Monograph 136, 1955, 5 p. To be republished in Inst. 
Elec. Engineers, Proc., Part C, size 27.5 cm. 

This table gives L,(x) to 5 or 6D for m = 0(1)10 and for x = 0(.1)5, together 
with modified second differences: 


5m? = 5 — 0.1846* + 0.03808268; 


it was computed on EDSAC 1 by direct summation of the series 
La(x) = X (—n)(—m +1) +++ (=m ter — 1) xt/ (rt 
r=0 


It is stated that the table was checked by hand differencing in the x-direction 
and by use of the recurrence relation 


(m + 1)Ln4. = (2m +1 — x)L, — nLa-1 


in the -direction. 

The table is stated to be correct to six places of decimals. However, for large 
n and x only 5 decimals are given for L,(x). Spot checking of the table, by calcu- 
lation from the explicit expression for L,(x) given in Admiralty Computing 
Service Report 82 [1 ] revealed non-trivial errors, e.g., L(4.3) is given as 0.08940 0 
whereas we find Ls = [3604.23764 321]/40320 = 0.08939 081--- giving a dis- 
crepancy of nine units in the sixth place. 

The author has discovered a systematic error and believes the listed values 
to be too large by one or more units in the sixth decimal for m = 5 and x > 4.1, 
n = 6 and x > 3.4, n = 7 and x > 2.8, m = 8 and x > 2.4, n = 9 and x > 2.1, 
and m = 10 and x > 1.9. 

The choice of a 6D multiplier for 6* and a 3D one for 5* seems curious, and 
can only add confusion to the transatlantic battle over .184 as against .18393. 
The multiplier .038 is suggested by L. J. Comrie [2]. 

Another unconventional point is that the differences are given, not in units 
of the last place of the tabulated values, but in full, and to fewer places than the 
main table. Thus we find: 


Lyo(5.0) = 1.75628 5,2 = —0.0500. 
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J. E. Wilkins and Nina Kropoff [3] give L,(x)/n! to 4D for m = 2(1)7, 
x = 0(.1)10(.2)20. 

There have been various tables of the corresponding orthogonal functions, 
e-*L, (x), or multiples thereof. For instance, F. G. Tricomi [4], who has studied 
the asymptotic behaviour of the L,,*(x), tabulated e-*L, (x) to4D for m = 1(1)10, 
x = .1(.1)1(.25)3(.5)6(1)14(2)34. This table, originally published in 1941 [5] has 
been reproduced in the 1948 edition of Jahnke-Emde [6]. A table has also been 
given by N. Wiener [7]. A table on punched cards, computed under the super- 
vision of C. Lanczos, is available at Numerical Analysis Research, University 
of California at Los Angeles, California. It gives e*L,(2x) to 10D for the fol- 
lowing ranges: m = 0(1)12, x = 0(1)12, x = 0(.1)5(5)10(1)20, m = 13(1)20, 
x = 0(.5)10(1)20. 

x 

1. Great Britain, H. M. Nauticat Atmanac Orrice, Zeros of Laguerre Polynomials and the 
— Christoffel Numbers, MTAC, v. 2, 1946-47, RMT 252, p. 31. 

2. L. J. Comrie, Chamber's Six-Figure Mathematical Tables, v. II, Natural Values, D. Van 
Nostrand C Co., Inc., ‘New York, 1949, p. 533. 


J. E Wixrns, Jr. & Nina ichosese, Table of Laguerre Functions, MTAC, v. 5, 1951, 
oy Le. p. 232. 
F. G. TRicoMt, “Sul comportamento asintotico dei polinomi di uerre,”” Ann. di Mat., 
ss. : v. 28, 1949, p. 263-289. [MTAC, v. 4, 1950, RMT 828, p. 216-217. 
5. Francesco Tricomi, ‘“Generalizzazione di una formula asintotica sui polinomi di Laguerre 
e sue applicazioni,”’ Accad. delle Scienze di Torino, Cl. d. sci. fis., mat., e nat., Atti, v. 76, 1941, 
P. 288-316. [MTAC, v. 2, 1947, RMT 385, p. 267. ] A revised and extended version of this table 
as just been published in the same Semanal (v. 90, 1955-56, p. 1-8); it will be reviewed in the 
next issue of this journal. 
6. Fritz Empe, Jahnke-Emde Tables of Higher Functions, Treated by Frits Emde. Fourth 
— ee with 177 figures. Leipzig, Teubner, 1948. [MTAC, v. 3, 1949, RMT 596, p. 
64-. 
7. NORBERT WIENER, Extrapolation, Interpolation, and Smoothing of yey Time Series 
with Engineering Applications, ass. Inst. Tech., Cambridge, Mass., and John Wiley & Sons, 
New York, 1949. [MTAC, v. 4, 1950, RMT 718, p. 24-25. ] 


1 z 
57(L ].—G. M. Mutter, Table of the Function Kj,(x) = - f u"Ko(u)du, Office 


of Technical Services, Department of Commerce, Washington 25, D. C., 1954, 
91 p., Price 70 cents. 


The tables include Kj, (x) = “ f u"Ko(u)du, n = 0(1)31, x = 0(.01)2(.02)5, 


7 to8D. Here Ko(u) is the modified Bessel function of the second kind of order zero. 

The tables were calculated on punched card equipment, and several significant 
zeros are printed after the last significant digit. They are accompanied by an 
introduction describing the means of calculation and drawing attention to some 
uses—particularly the calculation of the integral of the product of a function 
whose Taylor’s expansion is at hand and the function Ko(u). An example is 
worked out for illustration. 

No aids to interpolation are included. 

The tables were checked by computation which would have accumulated 
errors in the evaluation of Kjo(x), and the values for this table were checked 
against values already published by Bickley and Naylor [1]. No differencing or 
other such checks are mentioned. 

The entries are legible, but no spacing is used to improve legibility. 
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Related tables are listed in Fletcher, Miller, and Rosenhead, An Index of 
Mathematical Tables, sections 20.832 and 20.833. 


C. ‘BT. 


1. W. G. BicxLey & J. Naytor, “A short table of the functions Ki, (x) from m = 1 ton = 16,” 
Phil. Mag. (7), v. 20, 1935, p. 343-347. 


58[L ].—Joseru Kaye, ‘‘A table of the first repeated integrals of the error func- 
tion,” J. Math. and Phys., 1955, p. 119-125. 


The tables contain values of the repeated error integrals defined by the 
relations 


irerfex = f i" erfc dt, n =1,2,--- 


2 
t erfox = Pr ee 
Va 


and satisfying the recurrence relation 
(1) 2nyn = Yn-2 — 2XVn-1 


where y, = 7” erfc x. The values are given form = —1(1)11 forx = 0(.05)1(.1)2.8 
or 3.0 to 6D or 9D together with second central differences. The tables were 
computed with the aid of the recurrence relation (1) applied in the order of in- 
creasing m starting with the values i and 7° available to fifteen places [1] in 
the National Bureau of Standards tables. This extreme accuracy was necessary 
to obtain six significant figure accuracy in 7! because of the way in which the 
recurrence formula was employed. It is of interest to note that the functions 
might have been generated in a manner similar to that used [2] in the British 
Association tables. Thus if values of y, are desired for n = 0, 1, ---, p we rewrite 
(1) in the form 


(2) Yn-2 _ 2ngn + 2x9 n—1 


and start for nm > », with an arbitrary choice for 7, and J,-1, say Jn = 0, Jn-1 = 1 
and generate a sequence j, for k = n — 2, n — 3, ---,0, —1. If m is chosen 
sufficiently large then we obtain y, = af, where a is a constant, for k = —1, 0, 
1, ---,q where g < m. The constant may be simply determined from the relation 


2 
y.1 = —_— = af_1. Then increasing » (by 5, say) we again generate the 


sequence §,, comparing the computed values of y,., k = 0(1)p. 

This procedure can be repeated until the desired agreement is obtained. The 
technique is particularly suited to high-speed computers and has the advantage 
of producing all the desired function values without any loss of significant figures. 

Related tables are listed in Fletcher, Miller, and Rosenhead, An Index of 
Mathematical Tables, section 155. 


MILTON ABRAMOWITZ 
National Bureau of Standards 


Washington, D. C. 
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NBS Applied Mathematics Series No. 41, Tables of the Error Function and its Derivative, 
v.s Gov. Printing Office, Washington, D. C., 1954, xi + 302 p. Price $3.25. 
2. BRITISH ASSOCIATION MATHEMATICAL TABLES, v. X, Bessel Functions, Part 11. Functions 
of Positive Integer Order, Cambridge, 1952. 


59[P, S, T, Z].—Artuur Rose, R. Curtis Jonnson, Ricuarp L. Heiny, & 
THEODORE WILLIAMS, “Computers, Mathematics, Statistics, and Automa- 
tion,” Industrial and Engineering Chemistry, v. 48, 1956, p. 622-632. 


This paper is included in the fourth annual review of fundamentals of chemical 
engineering published each spring in Industrial and Engineering Chemistry. The 
authors attempt to give an exhaustive review of the year’s literature pertaining 
to applications of computers, mathematics, and statistics in engineering chem- 
istry. In addition, they discuss the current trends toward automation in chemical 
engineering in a way which provides a good background for chemical engineers 
interested in the subject. 

The present paper and its three predecessors in the series, [1], [2], [3], 
constitute a serious attempt to make available to chemists and other scientists 
and engineers the current status of applications of computers, mathematics, and 
statistics in the chemical industry. 

Cc. Be. 


1. ARTHUR Rose, Joan A. Scuitk, & R. Curtis Jounson, “Computers, Statistics, and Mathe- 
matics,” Industrial and Engineering Chemistry, v. 45, 1953, p. 933-939. 

2. ARTHUR Rose, RicHARD L. HeEtny, R. Curtis JoHNson, & Joan A. Scuiik, “Computers, 
Statistics, and Mathematics,” Industrial and Engineering Chemistry, v. 46, 1954, p. 916-922. 

3. ArTHUR Rose, R. Curtis Jonnson, & RicnHarp L. Hetny, “Computers, Statistics, and 
Mathematics,” Industrial and Engineering Chemistry, v. 47, 1955, p. 626-632. 


60[S].—NBS Applied Mathematics Series, No. 10, Tables for Conversion of 
X-ray Diffraction Angles to Interplanar Spacing, U. S. Govt. Printing Office, 
Washington, D. C., 1950, iii + 159 p., 27 cm. Price $1.75. 


The modern general purpose electronic computers, as they become more 
widely distributed, will, in time, render the use of tables largely obsolete. In 
reviewing the appearance of any new book of tables its proposed use should 
therefore be subjected to scrutiny to determine whether the existence of the 
tables is of lasting value or not. The interplanar spacings in the table being re- 
viewed are required in their own right, as, for instance, in chemical identification. 
For this purpose they are compared with spacings given in the ASTM index, and 
are not, in general, used as part of further complex calculations. It seems probable, 
therefore, that the tables will continue to be useful to the scientist for some 
time to come. 

The number of significant figures given is more than adequate for identifica- 
tion purposes. Indeed such accuracy for small values of the diffraction angle when 
the a-doublet is unresolved is unjustified in practice. The observed angle does not 
correspond to the X-ray wave-length Axa, used in making the tables, but to a 
weighted mean of Axa, and Axa,- For the high angle diffraction spectra used to 
find unit cell dimensions, the accuracy is justified, but the tables are incomplete 
since it is desirable to use the Kaz spectra as well as the Kai. 

The values of the interplanar spacings have been calculated for the six most 
commonly used wave-lengths in X-ray crystallography. It seems an unnecessary 
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luxury to repeat two of these six tables in terms of twice the angle and even if 
this is condoned, the choice of FeKa for one of the two preferentially treated 
wave-lengths can be queried. 

Nevertheless it is felt that the tables do satisfy adequately a present and 
continuing need, especially as such tables will not be given in the forthcoming 
volumes of the International Tables for X-ray Crystallography. 

A. S. DouGias 


A. M. B. DovuGtLas 
Cambridge 


England 


61[S].—G. E. Brown & D. F. Mayers, “The coherent scattering of y-rays by 
K electrons in heavy atoms,”’ R. Soc. London, Proc., v. 234A, 1956, p. 387-389. 


The matrix elements describing the coherent scattering of y-rays by electrons 
may be written as series involving associated Legendre polynomials. The coeffi- 
cients of the polynomials entering these series were evaluated to three decimal 
places on EDSAC, and it is shown that to this accuracy the coefficient of P? 
is negligible. 

i, Ti 


62[Z].—James R. NeEwMan, What is Science?, Edited and with an introduction 
by James R. Newman. Simon & Schuster, New York, 1955, viii + 493 p. 
Price $4.95. 


This is a collection of essays by twelve scientists in which they explain their 
fields to the layman; each is preceded by a biographical sketch and there is a 
bibliography and index. Although this review deals mainly with the contribution 
of Dr. J. Bronowski, entitled ‘‘Science as foresight,” it is appropriate to call atten- 
tion to two other articles, for their authors have made notable contributions to 
the area of mathematics with which this journal is concerned. The late Sir 
Edmund Whittaker, who discusses Mathematics and Logic, recalled that shortly 
after his appointment to the chair of mathematics at Edinburgh, he instituted the 
first University Mathematical Laboratory. Arising out of courses given there came 
his book with G. Robinson, The Calculus of Observations, which has been a 
basic text in classical numerical analysis. Dr. E. U. Condon, who writes on 
Physics, was responsible for the establishment of organizations of mathematicians 
and automatic computing machines devoted to modern numerical analysis. 

Dr. Bronowski, who is now Director of the Central Research Establishment 
of the (British) National Coal Board, is distinguished for contributions to such 
varied fields as algebraic geometry, operations research, literary criticism, as 
well as philosophy of science; he is also well known to the listening public in 
Europe for his talks and dramas. 

In his essay Bronowski discusses the progress of science as a development of 
tools which are extensions of the brain, rather than the hand. He begins with a 
discussion of deductive machines, underlining their speed and logical character 
and discusses the programming of two simple problems: nm? = 1+3+5+-::: 
+ (2m — 1), Xn41 = Xn(2 — ax,). The next section deals with adaptive machines. 
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He contrasts the case of machines programmed to play games using unlimited 
foresight (exhaustive search of all developments, or using an analytic theory) 
with those using limited foresight (exploring a short sequence of moves) and then 
accumulating experience. There is then a discussion of the machine of Ross Ashby. 

There follows a discussion on the theory of games exemplified by two-fingered 
Morra, and this leads to the ideas of “‘chance”’ and “random” and thence to those 
of “system” and “order.” The concepts of information and redundancy are then 
discussed. 

The final section deals with the logic of science. Bronowski suggests that while 
the more mechanical aspects of scientific method are comparable with the ac- 
tivities of deductive and adaptive machines, the characteristically human part 
is that of foresight based on insight, which is developed by a discovery of patterns. 
He then develops the idea of scientific theories as decoding: the aim of science 
is to present what we know in the most coherent fashion. 

It has been pleasing to find essays from newcomers as well as practiced ex- 
positors. The book will be specially appreciated by those of us who have no gift 
for the popular exposition of our field, and also by those of us who are associating 
ourselves and our machines with scientists who work in fields unfamiliar to us. 


he 


63[Z].—D. G. Prinz, “Programming the Transportation Problem for Digital 
Computers,”’ Technical Note CS65, Ferranti Ltd., Moston, Manchester, 
England, 1955, 32 p., mimeographed. 

The author describes the simplex method for solving linear programming 
problems, especially that version in which one transforms the inverse of the basis 
rather than the entire tableau, and then shows how this method may be coded 
to solve the transportation problem. He closes with some suggestions on the 
judicious use of ternary arithmetic for storing the inverse of the basis (which, for 
the transportation problem, has all of its entries 0, 1 or —1). 


A. J. HoFFMAN 
National Bureau of Standards 
Washington, D. C. 


64[Z ].—Harvey Coun, Code Library for ‘Stability Configurations of Electrons 
on a Sphere,” deposited in Library of Automatic Computer Codes. 15 tabular 
sheets, 37 cm. IBM Model 701 computer. 


These are coded machine commands for the calculations reported in [1]. 

Any one wishing to repeat the computation on the 701 would probably profit 
from one portion of the code (and probably only that one), the sub-routine for 
finding force vectors (No. 587 to 632) sheets 5, 6, and 7, and (No. 760 to 773) 
sheet 9. 


Harvey CoHuN 
Department of Mathematics 
Wayne University 
Detroit 1, Michigan 


1. Harvey Conn, “Stability configurations of electrons on a sphere”’ [p. 117 this issue]. 
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TABLE ERRATA 
In this issue reference has been made to errata in Review 56. 


248.—Haro tp T. Davis, Tables of the Higher Mathematical Functions, v. 1, The 
Principia Press, Inc., Bloomington, Indiana, 1933. 


The following errata have been found. 


The Gamma Function, Table 1, p. 201, x = 1.0255: 
For 0.98590 26815 read 0.98590 94917. 

The Gamma Function, Table 4, p. 250, Log I’ (22.7): 
For 20.5459 7344 7877 read 20.6459 7344 7877. 


A systematic error in Table 4 was discovered some time ago for Log I(x), 
x = 69.9 + n, where » = 0(1)31. The corrected values are as follows: 


p. 253 for Log I'(69.9) read 98.0491 3939 9253 p. 255 for 

p. 253 for Log I'(70.9) read 99.8936 1657 4999 Log '(86.9) read 130.1906 2460 9697 
p. 254 for Log '(71.9) read 101.7442 6281 0182 Log I'(87.9) read 132.1296 4438 6145 
p. 254 for Log 739) read 103.6009 9170 0565 Log '(88.9) read 134.0736 3326 1219 
p. 254 for Log I'(73.9) read 105.4637 1922 8883 Log I'(89.9) read 136.0225 3502 2189 
p. 254 for Log I'(74.9) read 107.3323 6366 7278 Log F'(90.9) read 137.9762 9471 3923 
p. 254 for Log 1'(75.9) read 109.2068 4548 4977 Log P'(91.9) read 139.9348 5859 7145 
p. 254 for Log ['(76.9) read 111.0870 8726 0872 Log P'(92.9) read 141.8981 7410 8531 
p. 254 for Log rye read 112.9730 1360 0674 Log P'(93.9) read 143.8661 8982 2524 
p. 254 for Log I'(78.9) read 114.8645 5105 8346 Log 1'(94.9) read 145.8388 5541 4790 
p. 254 for Log I'(79.9) read 116.7616 2806 1556 Log P'(95.9) read 147.8161 2162 7218 
p. 254 for Log a read 118.6641 7484 0870 Log O79 read 149.7979 4023 4388 
p. 254 for Log I'(81.9) read 120.5721 2336 2482 Log P'(97.9) read 151.7842 6401 1439 
p. 254 for Log 3 read 122.4854 0726 4243 Log P'(98.9) read 153.7750 4670 3242 
p. 254 for Log I'(83.9) read 124.4039 6179 4793 Log T1002 read 155.7702 4299 4839 
p. 254 for Log I'(84.9) read 126.3277 2375 5621 Log '(100.9) read 157.7698 0848 3065 
p. 254 for Log P'(85.9) read 128.2566 3144 5865 


The first of these errata was reported by Professor Charles A. Hutchinson, 
University of Colorado, Boulder, Colorado. The others were furnished by Pro- 
fessor H. T. Davis. 
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NOTES 


Alan Mathison Turing 
1912-1954 


Dr. A. M. Turing, who played a decisive part in various phases of the develop- 
ment and exploitation of automatic computing machines, died on 7 June 1954; 
he was born in London on 23 June 1912. A detailed account of his life and work, 
which we have used in the preparation of this note, has been prepared by M. H. 
A. Newman [1]. 

It was about 1935 that Turing, first at Cambridge and then at Princeton, 
began studies in mathematical logic which led him to introduce the concepts of 
“computable numbers” and what are now known as “Turing machines.”’ It was 
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not until after the war that he became active in the realization of machines of 
this type. During 1945-7, when Turing was at Teddington, laying the foundations 
for the Automatic Computing Engine of the National Physical Laboratory, we 
were living close by. We recall visits from Turing, with long discussions on various 
mathematical topics; we always found him an active listener, willing and able to 
contribute ideas on many fields, some far away from his current interest. In 
addition, from our window, we often saw him passing by on his training runs. 
From 1948 until his death he was at the University of Manchester participating 
in the mathematics and machine development program there. For his services to 
his country in World War II, he was awarded the O.B.E.; he was elected to the 
Fellowship of the Royal Society in 1951. 

The classical paper, ‘On computable numbers, with an application to-the 
Entscheidungsproblem,’’ London Math. Soc., Proc., s. 2, v. 42, 1936, p. 230-265; 
s. 2, v. 43, p. 544-546, was followed by one indicating the equivalence of his 
concept of “computable” with those of “general-recursive” (Gédel-Kleene) and 
“\-definable’”’ (Church). Among his other papers in logic was one on the insolu- 
bility of the word-problem for semi-groups with cancellation. In this connection 
we point out his excellent expository article : ‘Solvable and unsolvable problems,” 
Science News, v. 31, 1954, p. 7-23. 

Turing was responsible for the initial version of ‘Programmers’ Handbook of 
the Manchester machine,” 1950. His published work in numerical analysis was in 
two directions. In the first, “‘Rounding-off errors in matrix processes,”’ Quart. Jour. 
Mech. and Appl. Math., v. 1, 1948, p. 287-308, he introduced the concept of 
“condition-number,” which has played a large part in later studies. In the second 
he was concerned with the Riemann zeta-function. His pre-war theoretical in- 
vestigations were described in ‘‘A method for the calculation of the zeta-function,” 
London Math. Soc., Proc., s. 2, v. 48, 1943, p. 180-197. Some computing was 
done in 1950 on the Manchester machine; in the resulting paper, ‘Some calcula- 
tions of the Riemann zeta-function,”’ London Math. Soc., Proc., s. 3, v. 3, 1953, 
p. 99-117, we find an unusually revealing story of work in progress and valuable 
remarks on the idea of rigorous computation. 

Two other very readable articles are, ‘Computing machinery and intelligence,” 
Mind, v. 59, 1950, p. 433-460, and ‘“‘Digital computers applied to games: chess,” 
Faster than Thought, ed., B. V. Bowden, Pitman, London, 1953, p. 288-295. 

His latest work was in the domain of applied mathematics, and was concerned 
with the problem of morphogenesis, that is, of the growth and form of living 
material. He developed a chemical theory for this phenomenon and found it 
satisfactory, when tested on the Manchester computer, in certain botanical situa- 
tions; this is described in ‘“The chemical basis of morphogenesis,’ Roy. Soc., 
Phil. Trans., Sec. B, v. 237, 1952, p. 37-72. 

OLGA TAUSSKY 


oo? A 
National Bureau of Standards 


Washington, D. C. 
1. Biographical Memoirs of Fellows of the Royal Society, v. 1, 1955, p. 253-263. 











182 CORRIGENDA 
CORRIGENDA 


N. G. W. H. BEEcEr, Table of the Least Factor of the Numbers that are not 
Divisible by 2, 3, 5 of the Eleventh Million, MTAC, v. 10, 1956, Review 5, p. 36. 

Lines 2 and 3, for 429 tables, 34 pages of form, 333 cm., deposited in the UMT 
FILE, read 462 p. including 429 p. of tables on 33 fascicles, each fascicle having 
the same introductory page except for indication of serial number and limits, 
333 cm., deposited in the UMT Fire. 


GIUSEPPE PALAMA & L. Po ettt, “‘Tavola dei numeri primi dell’intervallo 
12 012 000-12 072 060,” MTAC, v. 10, 1956, Table Errata 247, p. 54. 
Line 5, for Division read Divisor. 
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